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Appendix:  Listing  and  subroutines  details 

One  sample  program  (with  its  output)  and  18  subroutines  are  listed  in  this  appendix. 


•  sample  program  and  its  output 

•  subroutine  cschur 

•  subroutine  cisol 

•  subroutine  cmhu 

•  subroutine  chqr 

•  subroutine  corder 

•  subroutine  cexphy 

•  subroutine  cexpri 

•  subroutine  cindex  (in  cexpri) 

•  subroutine  cmswap  (in  cexpri) 

•  subroutine  cexpnt 

•  subroutine  cddexp 

•  subroutine  cfmulv 

•  complex  function  cdotc  (Unpack  BLAS,  in  blas_c) 

•  complex  function  cdotu  (Unpack  BLAS,  in  blas_c) 

•  subroutine  caxpy  (Unpack  BLAS,  in  blas_c) 

•  subroutine  ccopy  (Unpack  BLAS,  in  blas_c) 

•  subroutine  csscal  (Unpack  BLAS,  in  blas_c) 

•  subroutine  cscal  (Unpack  BLAS,  in  blas_c) 
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Sample  program 


THIS  PROGRAM  SHOWS  THE  USAGE  OP  SUBROUTINES  FOR  COMPUTING  THE  MATRIX 
EXPONENT  I AL . 

PARAMETER  (M=6) 

COMPLEX  P(M,M)  ,  B(M,M)  ,E(M,M)  ,W(5*M)  ,TAU,V(M,M) 

REAL  OVFT 

INTEGER  M,  N,  I  , J , IP(M)  , IDX(M)  .  I  ERR , NP , NB ,NE , NV 
C  A  NUMBER  NEAR  THE  OVERFLOW  THRESHOLD  ON  A  VAX  IN  SINGLE  PRECISION 
OVFT=l  0E3  6 

C  ROV  DIMENSION  OF  THE  MATRICES 
NB=M 
NP=M 
NE=M 
NV=M 

C  READ  IN  THE  DIMENSION,  TAU .  AND  THE  WHOLE  MATRIX  (BY  ROWS) 

READ*,  N,  TAU.  ( ( B(  I  . J )  , J=1  . N)  ,  1=1  , N) 

PRINT*,  ’ INPUT  MATRIX’  ,  "TAU="  ,  TAU 
CALL  MATPT(M.N.B) 

C  SET  V=l 

DO  2  0  J=1 , N 
DO  10  1=1 ,N 
10  V ( I , J )=CMPLX( 0 ) 

20  V ( J , J )=CMPLX ( 1 ) 

C  COMPUTE  THE  MATRIX  EXPONENTIAL 

CALL  CSCHUR  ( 8 , P ,  W,  N , NB , NP ,  I  ERR ) 

CALL  CORDER ( B , P , TAU , W, N , NB , NP ) 

PRINT*,  ’  ’ 

PRINT*,  *  SCHUR  FORM  S* 

CALL  MATPT ( M , N , B ) 

PRINT*,  '  ’ 

PRINT*,  "UNITARY  MATRIX  P  * 

CALL  MATPT (M,N , P ) 

CALL  CEXPHY ( B , E ,  W,  IP,  I DX , TAU , OVFT ,  I  ERR , N , NB , NE ) 

I F (  I  ERR  EQ  -2)  THEN 

PRINT*,  "THE  EXPONENTIAL  OF  SOME  EIGENVALUE  OVERTLOVS* 

STOP 
END  IF 

CALL  CFMUL V( E , P , V ,  W,  N , NE , NP , NV , N ) 

PRINT* ,  ”  " 

PRINT*,  " EXP ( TAU • B ) " 

CALL  MATPT (M, N , V ) 

END 

SUBROUTINE  MATPT (M.N, X) 

COMPLEX  X(M,N) 

DO  10  1=1  ,N 

WR I TE( 6 ,  10  0)  ( X (  I  , J )  , J=1  ,N) 

10  CONTINUE 

100  FORMAT) IX, 3 ( " ( " , 2G1 2  3,*)")) 

RETURN 

END 
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We  ran  the  program  on  a  V ax-780  and  obtained  the  following  results.  For  simplicity,  only 
three  digits  are  displayed,  and  the  output  has  been  reformatted  so  that  it  is  easy  to  read. 
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CSCHUR 


1.  PURPOSE 

.  The  Fortran  77  subroutine  CSCHUR  reduces  a  complex  general  matrix  B  to  a  complex 
upper  Triangular  matrix  S  using  unitary  similarity  transformations. 


2.  USAGE 

(A).  Calling  Sequence. 

SUBROUTINE  CSCHUR(B,P,w,n,nb,np,ierr) 

Parameters  : 

B  is  a  complex  two-dimensional  variable  with  row  dimension  nb  and 
column  dimension  at  least  n.  On  input,  B  contains  the  complex  matrix 
of  order  n  to  be  reduced  to  triangular  form.  On  output,  B  is  overwritten 
by  its  Schur  form  S  (upper  triangular). 

P  is  a  complex  two-dimensional  output  variable  with  row  dimension  np 
and  column  dimension  at  least  n.  On  output,  P  contains  the  unitary 
matrix  that  transform  B  to  S  :  S=P*  B  P. 

w  is  a  complex  one-dimensional  variable  with  dimension  2n.  Work  space. 

n  is  an  input  integer  equal  to  the  order  of  the  matrix  B. 

nb,np  are  input  integers  equal  to  the  row  dimension  of  B  and  P  respectively. 

ierr  is  an  integer  variable  indicating  the  error  condition.  If  the  number  of 
QR-iterations  for  some  eigenvalue  reaches  30  without  converging  (at 
that  point  the  subroutine  CHQR  terminates),  then  ierr  is  set  equal  to 
-1;  otherwise  ierr=0. 


Subprograms: 


CISOL,  CMHU,  CHQR. 


3.  DISCUSSION  OF  METHOD  AND  ALGORITHM 

We  first  call  CISOL  to  isolate  the  eigenvalues  of  B.  Then  we  employ  Householder 
transformations  to  reduce  B  to  upper  Hessenberg  form  H  (subroutine  CMHU).  Finally 
we  apply  QR-algorithm  to  transform  H  to  upper  triangular  form  (subroutine  CHQR). 
For  details  of  the  algorithm,  see  the  description  of  subroutines  CISOL,  CMHU  and 
CHQR. 

There  are  approximately  3  executable  Fortran  statements. 
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5.  PROGRAM  LISTING. 


SUBROUTINE  CSCHUR  ( B , P ,  W,  N , NB , NP ,  l ERR ) 

COMPLEX  B(NB,N) ,P(NP,N) ,  W(  2  •  N ) 

INTEGER  N.NB.NP. I  ERR 

CSHCUR  FIRST  ISOLATES  THE  EIGENVALUES  OF  B  WHENEVER  POSSIBLE.  THEN  USES 
HOUSEHOLDER  TRANSFORMATION  TO  REDUCE  B  TO  UPPER  HESSENBERG  H,  FINALLY 
USES  QR  TO  TRANSFORM  H  TO  ITS  SCHUR  FORM  S  THAT  WILL  REPLACE  B  THE 
TRANSFORMATIONS  ARE  ACCUMULATED  IN  P,  IE  ,  S=PH»B*P,  HERE  PH  STANDS 
FOR  THE  CCMPLEX  CONJUGATE  TRANSPOSE  OF  P  IF  QR  FAILS  TO  CONVERGE.  WE 
RETURN  I ERR=-  1 ,  OTHERWISE  I ERR=0 

CODE  IN  F  7  7  BY  K  C.  NG ,  UC  BERKELEY,  REVISED  ON  5/22/85 
REQUIRED  SUBROUTINES  :  Cl  SOL,  CMHU ,  CHQR 

GLOSSARY 

B  INPUT  A  COMPLEX  MATRIX,  OUTPUT  TRIANGULAR  MATRIX 

P  A  UNITARY  MATRIX 

W  CCMPLEX  ARRAY,  WORK  SPACE 

N  DIMENSION  OF  MATRIX  H 

NB , NP  LEADING  DIMENSION  OF  ARRAys  B  AND  P  RESPECTIVELY 
I  ERR  OUTPUT  (INTEGER)  ERROR  STATUS  IF  QR  DOESN’T  CONVERGE  IN  30 
ITERATION,  THE  PROGRAM  QUITS  AND  RETURNS  I ERIfc=- 1 

INTERNAL  VARIABLES 

INTEGER  LCW,  IGH 

ISOLATE  THE  EIGENVALUES  IN  B 

CALL  CISOL(B,P, LCW, IGH, N.NB.NP) 

TRANSFORM  B  TO  UPPER  HESSENBERG  FORM 

CALL  CMHU ( B, P, LOW, IGH, N,NB,NP) 

TRANSFORM  THE  UPPER  HESSENBERG  MATRIX  TO  TRIANGULAR  FORM  AND  ACCUMULATE 

CALL  CHQR ( B , P , LCW,  JGH,W,W{N+l)  ,N,NB,NP,  I  ERR ) 

RETURN 


cschur-cisol 


CISOL 

(Subprogram  of  CSCHUR) 


1.  PURPOSE 

The  Fortran  77  subroutine  CISOL  isolates  (using  permutations)  the  eigenvalues  of  B 
whenever  possible. 


2.  USAGE 

(A).  Calling  Sequence. 

SUBROUTINE  CISOL(B,P,low,igh,n,nb,np) 

Parameters  : 

B  is  a  two-dimensional  complex  variable  with  row  dimension  nb  and 
column  dimension  at  least  n.  On  input,  B  contains  a  complex  matrix  of 
order  n.  On  output,  B  contains  the  permuted  matrix. 

P  is  an  output  two-dimensional  complex  variable  with  row  dimension  nb 
and  column  dimension  at  least  n.  P  contains  the  permutations. 

low.igh  are  integer  output  variables  indicating  the  boundary  indices  for  the  per¬ 
muted  matrix. 

n  is  an  input  integer  equal  to  the  column  dimension  of  the  matrix  B. 

nb,np  are  input  integers  equal  to  the  row  dimension  of  B  and  P  respectively. 


3.  DISCUSSION  OF  METHOD  AND  ALGORITHM 

This  is  a  F77  variation  of  the  first  half  (isolating  the  eigenvalues  of  B  whenever  possi¬ 
ble)  of  cbal,  an  EISPACK  Fortran  IV  subroutine.  The  permutations  are  recorded  in  a 
matrix  P.  For  details  of  the  algorithm,  see  EISPACK  [2)  or  Parlett  and  Reinscb  ( 1] . 

There  are  approximately  42  executable  Fortran  statements. 


4.  REFERENCES 

jl]  B.N.  Parlett  and  C.  Reinsch,  Balancing  a  Matrix  for  Calculation  of  Eigenvalues 
and  Eigenvectors,  Num.  Math.  13,  293-304  (1969).  (Reprinted  in  Handbook  for 
Automatic  Computation,  Vol.  II,  Linear  Algebra,  J.H.  Wilkinson  -  C,  Reinsch, 
Contribution  II/ 1 1,  315-326,  Springer-Verlag,  1971.) 

[2|  B.T.  Smith  &  et  al,  Matrix  Eigensystem  Routines  -  EISPACK  Guide,  Lecture 
Notes  in  Computer  Science,  Vol  6,  Springer-Verlag,  1974. 
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cschur-cisol 


5.  PROGRAM  LISTING. 


SUBROUTINE  C I  SOL ( B , P . LOW,  IGH.N.NB.NP) 

COMPLEX  B(NB.N) . P(NP.N) 

INTEGER  NB , N , LCW, I GH , NP 

C 1  SOL  ISOLATES  THE  EIGENVALUES  IN  THE  INPUT  COMPLEX  MATRIX  B.  IT  IS  A  F7 
VARIATION  OF  THE  FIRST  HALF  OF  EISPACK  SUBROUTINE  CBAL ,  WITH  THE  RESULTED 
PERMUTATION  RECORDED  IN  P  P  WILL  BE  SET  TO  AN  IDENTITY  MATRIX  INITIAL* 
CODE  IN  F77  BY  K  C  NG ,  UC  BERKELEY;  REVISED  ON  5/22/85. 

GLOSSARY 

B  A  COMPLEX  MATRIX 

P  ON  OUTPUT,  AN  UNITARY  MATRIX 

IGH.LCW  INTEGER  INDICES  INDICATING  BOUNDARIES 

N  DIMENSION  OF  MATRIX  B 

NB , NP  LEADING  DIMENSION  OF  ARRAYS  B  AND  P  RESPECTIVELY 


INTERNAL  VARIABLES 

INTEGER  I , J , K , L ,M, I EXC 

COMPLEX  X, ZERO 

DATA  ZERO/ ( 0E0 , 0E0 ) / 

SET  P=I 

DO  10  J=1  ,  N 
DO  5  1=1 ,N 
P ( I . J)=ZERO 
0  P(J,J)=1.0 


K=1 

L=N 

GOTO  100 

IN-LINE  PROCEDURE  FOR  ROW  AND  COLUMN  EXCHANGE 

0  CONTINUE 

IF(J.NE.M)  THEN 
DO  30  1=1 , L 
X=B( I , J ) 

B (  I  .  J  )=B(  I  ,M) 

B (  I  ,M)=X 

0  CONTINUE 

DO  <0  I=K,N 
X=B  ( J ,  I ) 

B(  J  ,  I  )=B(M,  1  ) 

B(M, I )=X 

0  CONTINUE 

DO  60  1=1 ,N 
X=P(  I  , J  ) 

P( I , J )=P( | ,M) 

P( I ,M)=X 

0  CONTINUE 

END  IF 

I F (  I EXC  EQ  2 )  GOTO  1 3  0 

SEARCH  FOR  ROWS  ISOLATING  AN  EIGENVALUE  AND  PUSH  THEM  DOWN 


100 


I  F ( L  EQ  1 )  GOTO  280 
L=L  -  1 

DO  12  0  J=L  .1.-1 


o  o  o 


cschur-cmhu 


CMHU 

(Subprogram  of  CSCHUR) 


l.  PURPOSE 

The  Fortran  77  subroutine  CMHU  reduces  a  complex  general  matrix  B  to  a  complex 
upper  Hessenberg  matrix  H  using  unitary  similarity  transformations.  The  Hessenberg 
matrix  H  is  needed  in  subroutine  CHQR  to  find  the  Schur  form  of  B. 


2.  USAGE 

(A).  Calling  Sequence. 

SUBROUTINE  CMHU(B,P,Iow,igh,n,nb,np) 

Parameters  . 

B  is  a  complex  two-dimensional  variable  with  row  dimension  nb  and 
column  dimension  at  least  n.  On  input,  B  contains  the  complex  matrix 
of  order  n  to  be  reduced  to  Hessenberg  form.  On  output,  B  is  overwrit¬ 
ten  by  the  upper  Hessenberg  matrix  H. 

P  is  a  complex  two-dimensional  variable  with  row  dimension  np  and 
column  dimension  at  least  n.  On  output,  P  contains  the  unitary  matrix 
that  transform  B  to  H  :  H— P*  B  P. 

low.igh  are  integer  input  variables  indicating  the  boundary  indices  for  the 
matrix.  See  CISOL  for  details. 

n  is  a  input  integer  equal  to  the  order  of  the  matrix  B. 

nb.np  are  input  integers  equal  to  the  row  dimension  of  B  and  P  respectively. 


S.  DISCUSSION  OF  METHOD  AND  ALGORITHM 

The  method  employs  Householder  transformations  to  reduce  B  to  upper  Hessenberg 
form.  A  general  description  (for  complex  B)  can  be  found  in  [1],  pp. 137-138.  For  real 
matrices,  EISPACK  [3]  (see  also  [4j )  has  a  subroutine  orthes  to  do  the  reduction.  The 
subroutine  CMHU  given  here  is  essentially  an  extension  of  orthes  to  complex  matrices. 
(EISPACK  has  another  subroutine  elmhes  to  reduce  a  general  complex  matrix  to 
Hessenberg  form,  but  the  transformations  are  non-unitary.) 

The  details  of  the  method  can  be  illustrated  by  a  6  by  6  example.  A  typical  stage  in 
the  reduction  would  be 
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where  the  elements  h  are  those  of  the  final  Hessenberg  matrix,  and  the  x’s  ,  t  and  f  are 
in  their  intermediate  states.  Let  T  and  F  denote  the  3  by  3  sub-matrices  of  B  consist¬ 
ing  of  small  letters  t  and  f  respectively.  Also  let  x  denote  the  vector  (xl,x2,x3)t,  as 
shown  in  the  above  figure.  The  first  step  of  the  Householder  transformation  is  to  find  a 
vector  u  to  form  the  Householder  matrix  R  such  that 


xl 

w 

R- 

x2 

= 

0 

x3 

0 

where  R  :=  I- 


uu 

~ 


d=u 


V2; 


then  the  sub-matrices  F  and  T  of  B  are  updated,  and  the  transformation  is  accumu¬ 
lated  in  P  as  indicated  : 

F  :=  R  F  R, 

T  :=  T  R, 


The  determination  of  w  and  the  vector  u  is  as  follows.  For  xl^O,  it  is  easy  to  show 
that 

W==±w^’  S;3=|xi|2+lx2l2+|x3|2. 

When  xl=0,  w  can  be  any  complex  number  with  magnitude  \/S.  In  CMHU,  we 
choose  w  to  be  (xl/|xl|)>/S  if  xl^O  and  VS  if  xl=0.  The  vector  u  is  just  equal  to 
(xl-w,x2,x3)‘.  However,  care  must  be  taken  in  computing  the  first  component 
U!=xl-w  when  xl  and  w  are  almost  equal  (cf.  [2],  p.91).  The  following  formula  is  used 
whenever  |w/xl|<2: 

.  -( |x2|z+  1  x3 1 2) 

U!  =  Xl-W  =  -1-1 - - 1 - 1— 

(xl+w) 


The  above  steps  are  repeated  on  further  columns  of  the  transformed  B  until  it  reaches 
the  Hessenberg  form. 

There  are  approximately  45  executable  Fortran  statements. 


4.  REFERENCES 

[l|  C.E.  Froberg,  Introduction  to  Numerical  Analysis,  2nd  ed..  Addison-  Wesley 
1969. 


(2|  B.N.  Parlett,  The  Symmetric  Eigenvalue  Problem,  Prentice-Hall,  Englewood 


cschur-cmhu 


Cliffs,  N.J.  1980. 

B.T.  Smith  &  et  al,  Matrix  Eigensystem  Routines  -  EISPACK  Guide,  Lecture 
Notes  in  Computer  Science,  Vol  6,  Springer-Verlag,  1974. 

J.H.  Wilkinson  &  C.  Reinsch,  Handbook  in  Automatic  Computation  Vol.  II, 
Linear  Algebra  in  Numerical  Analysis,  Springer-Verlag,  1971. 
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cschur-cmhu 


5.  PROGRAM  LISTING. 


SUBROUTINE  CMHU( B , P , LCW,  1 GH , N , NB , NP ) 

COMPLEX  B(NB,N) ,  P(NP,N) 

INTEGER  LCW,  IGH ,N,NB ,NP 

THIS  SUBROUTINE  REDUCES  A  CCMPLEX  MATRIX  B  TO  UPPER  HESSENBERG  FORM  BY 
UNITARY  TRANSFORMATIONS  TRANSFORMATIONS  ARE  ACCUMULATED  IN  P.  THIS 
SUBROUTINE  SHOULD  BE  PRECEEDED  BY  C I  SOL . 

CODE  IN  F77  BY  K.C.  NG .  UC  BERKELEY;  REVISED  ON  5/22/85. 

F77  GENERIC  FUNCTIONS  ABS  ,CONJG,  IMAG,  REAL,  SQRT 

F77  NON-GENERIC  FUNCTION:  CMPLX 

REQUIRED  BLAS  SUBROUTINES.  CAXPY ,  CCOPY ,  CSCAL ,  CSSCAL 
REQUIRED  BLAS  FUNCTION  CDOTC 

GLOSSARY: 

B  A  COMPLEX  MATRIX 

P  A  UNITARY  MATRIX 

I GH , LOW  INTEGER  INDICES  INDICATING  BOUNDARIES 

N  DIMENSION  OF  MATRIX  B,  P 

NB , NP  LEADING  DIMENSION  OF  ARRAyj  B  AND  P  RESPECTIVELY 


INTERNAL  VARIABLES 

CCMPLEX  X , Y , 2 , ZERO , CDOTC 
REAL  T , D , HALF , SCALE 
INTEGER  I , J , R , RP 1 

ZERO=0 
HALF=0 . 5 

IF  (IGH-LCWLT  2)  RETURN 
DO  2  00  R=LOW,  IGH- 2 

SCALE  COLUMN  AND  SET  UP  THE  HOUSEHObDER  VECTOR  IN  B(  ,R) 

RP 1=R+ 1 
SCALE=0 

DO  20  !=RP 1 , IGH 

0  SCALE=SCALE+ABS ( REAL ( B ( I ,R) ) )+ABS( I MAG ( B ( I , R ) ) ) 

I F( SCALE  EQ  REAL ( 2ERO) )  RETURN 
Y=0 

DO  SO  I=R+2 , IGH 
B( I ,R)=B( I ,R) /SCALE 
0  Y=Y+  B  (  I  ,  R  )  •  CON  JG  (  B  (  I  ,  R  )  ) 

Y=CMPLX( REAL ( Y ) ) 

IF(Y  EQ  ZERO)  GOTO  200 
X=e(RPl , R ) / SCALE 
I F ( X  EQ  ZERO)  THEN 
D=REAL ( Y ) 

Z=CMPLX( SQRT ( REAL (Y) ) ) 

B ( RP 1  , R )=-  Z 

ELSE 

T=SQRT ( REAL ( Y+X*CONJG( X ) ) ) /ABS(X) 

Z=X*CMPLX(T) 

I F ( T  LT  HALF)  THEN 
B( RP 1 ,R)=X- Z 

ELSE 

CANCELLATION  OCCURS,  USE  RATIONALIZED  FORMULA  FOR  X-Z 


c  s  c  fa  u  r - cmh u 


B( RP 1 , R)=-Y/CONJG(X+Z) 

END  IF 

D=REAL ( Y+B(  RP  1 , R ) *CONJG( B( RP1 , R ) ) ) ‘HALF 
END  IF 
C 

C  PERFORM  B-((B»U)/D)*UT,  P  -  ( ( P • U ) / D ) • UT .  UT  IS  THE  CONJUGATE  TRANSPOSE  OB 
C 

DO  9  0  1=1  ,N 

X-=0 
Y=0 

DO  70  J=RPl , IGH 

X=X+B(  1  . J  )  • B (  J . R ) 

Y=Y+P(  I  , J  )  • B ( J . R ) 

70  CONTINUE 

X=X/D 
Y=Y/D 

DO  80  J =RP! , IGH 

B (  I  ,  J )  =B  ( I  , J)-CONJG(B( J , R )  ) *X 
P (  I  , J )=P  ( I  , J ) - CONJG ( B (  J , R  )  ) • Y 
80  CONTINUE 

90  CONTINUE 

C 

C  PERFORM  B-U* ( ( UH*B) /D) 

C 

DO  120  J=RP1 , N 

X=CDOTC( IGH-RPl+1 ,6(RPl ,R) , 1 , B ( RP 1 , J) , 1 ) 

X=X/D 

CALL  CAXPY ( IGH-RPl+1 . -X , B( RPl ,R) , 1 , B( RPl , J ) , 1 ) 

1 2  0  CONT I NUE 

B( RPl , R )—Z  • SCALE 
C 

C  CLEAR  THE  ELEMENTS  BELOW  THE  SUB-DIAGONAL 
C 

DO  130  I=ft+  2 ,  IGH 
130  B( I ,R)=0 

200  CONTINUE 
RETURN 
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CHQR 

(Subprogram  of  CSCHUR) 


1.  PURPOSE 

The  Fortran  77  subroutine  CHQR  reduces  a  complex  upper  Hessenberg  matrix  H  to  a 
complex  upper  triangular  matrix  S  (Schur  form)  using  the  QR-algorithm.  The  unitary 
transformations  are  accumulated  in  P. 


2.  USAGE 

(A).  Calling  Sequence. 

SUBROUTINE  CHQRIH.P.low.igh.w.v.n.nh.np.ierr) 

Parameters: 

H  is  a  complex  two-dimeusional  variable  with  row  dimension  nh  and 
column  dimension  at  least  n.  On  input,  it  contains  the  complex  upper 
Hessenberg  matrix  H.  On  output,  it  contains  the  triangular  matrix  S, 
the  Schur  form  of  H. 

P  is  a  complex  two-dimensional  variable  with  row  dimension  np  and 
column  dimension  at  least  n.  On  input,  P  contains  an  unitary  matrix 
that  changed  B  to  H  :  H=»P*  B  P  (obtained  from  subroutine  CMHU). 
On  output,  P  contains  the  unitary  matrix  that  transforms  B  to  its 
Schur  form  S  :  S=P*  B  P. 

low.igh  are  integer  input  variables  indicating  the  boundary  indices  for  the 
matrix.  See  subroutine  CISOL  for  details. 

w,v  are  two  complex  one-dimensional  variables  of  order  n.  Work  space. 

n  is  an  input  integer  equal  to  the  order  of  the  matrix  H. 

nh.np  are  input  integers  equal  to  the  row  dimension  of  H  and  P  respectively. 

ierr  is  an  integer  variable  indicating  the  error  condition.  If  the  number  of 
QR-iterations  for  some  eigenvalue  reaches  30  without  getting  conver¬ 
gent  (the  subroutine  terminates  if  that  happens),  then  ierr  is  set  equal 
to  -1;  otherwise  ierr  =  0. 


S.  DISCUSSION  OF  METHOD  AND  ALGORITHM 

CHQR  is  essentially  a  complex  version  of  the  EISPACK  (3,4j  subroutine  hqr  (or  A^r  2) 
with  some  modifications.  The  method  is  to  apply  the  QR-aigorithm  (lj  to  the  Hessen¬ 
berg  matrix  H  to  reduce  it  to  triangular  form.  Since  the  QR-algorithm  is  a  popular 
topic  in  many  numerical  analysis  texts  (eg.  [2] ,(4) ),  we  will  not  repeat  the  general 
description  here.  Rather,  we  discuss  certain  technical  details  which  are  different  from 
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EISPACK’s  hqr  in  the  implementation  of  the  QR-algorithm.  There  are  four  major 
differences. 

(1) .  First  of  all,  complex  arithmetic  is  used  throughout  CHQR;  therefore,  an  orthogo¬ 
nal  matrix  becomes  an  unitary  matrix  as  described  in  [l|  and  double  shifting  for  conju¬ 
gate  eigenvalues  is  no  longer  needed. 

(2) .  The  criterion  for  testing  a  negligible  sub-diagonal  is  as  follows.  Let  |x|  stand  for 
|Re(x)|+ |Im(x)|,  and  let  t  denote  the  machine  precision.  Then  H(i,i-1)  is  negligible  if 
|H(i,i-l)| 

satisfies  (a)  and  (bl),  or  (a)  and  (b2)  if  (bl)  fails: 

(a).  |H(i,i-l)|  <  e  min{|H(i,i)|,|H(i-l,i-l)|}  ; 

(bl).  4|H(i,i-l)»H(i-l,i)l  <  t(|H(i,i)-H(i-l,i-l)  | 2); 

(b2).  2v'|H(i,i)*H(i-l,i-l)|  <  £-|H(i,i)+  H(i-l,i-l)|  . 

These  criteria  are  obtained  by  forcing  the  eigenvalues  of  the  following  2x2  sub-matrix 
to  agree  in  working  precision  with  the  eigenvalues  of  the  same  matrix  with  H(i,i-1) 
replaced  by  zero. 

f  H(i-l.i-l)  H(i-l.i)  | 

(  H(i.i-l)  H(i,i)  I' 


(3) .  CHQR  doesn't  look  for  two  consecutive  small  sub-diagonal  elements. 

(4) .  The  transformations  are  accumulated  in  P. 

There  are  approximately  70  executable  Fortran  statements. 
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5.  PROGRAM  LISTING. 


SUBROUTINE  CHQR  ( H , P , LOW,  IGH ,W,  V , N, NH , NP ,  I  ERR ) 

COMPLEX  II ( NH  , N )  ,  P  ( NP  ,N)  ,W(N)  ,V(N) 

INTEGER  LCW,  IGH.N.NH.NP,  I  ERR 

THIS  SUBROUTINE  USES  QR- ALGOR  ITHM  TO  TRANSFORM  A  COMPLEX  UPPER  HESSENBEB 
(MATRIX  H  TO  ITS  SCHUR  FORM.  THE  TRANSFORMATIONS  ARE  ACCUMULATED  IN  P.  THE 
SUBROUTINE  FAILS  (RETURN  I  ERR  =  -1  )  IF  ANY  EIGENVALUE  TAKES  MORE  THAN  8 
ITERATIONS . 

CODE  IN  F77  BY  K.C.  NG ,  UC  BERKELEY;  REVISED  ON  5/22/85. 

F77  GENERIC  FUNCTIONS  ABS  ,  CONJG ,  IMAG,  MAX,  REAL,  SQRT 

F77  NON-GENERIC  FUNCTION:  CMPLX 

GLOSSARY : 

H  INPUT:  A  COMPLEX  HESSENBERG  MATRIX;  OUTPUT:  TRIANGULAR  MATRIX 

P  AN  UNITARY  MATRIX 

W, V  COMPLEX  ARRAIES,  WORK  SPACE 

IGH.LCW  INTEGER  INDICES  INDICATING  BOUNDARIES 

N  D (MENS  ION  OF  MATRIX  H 

NH , NP  LEADING  DIMENSION  OF  ARRAIES  H  AND  P  RESPECTIVELY 
I  ERR  OUTPUT  (INTEGER)  ERROR  STATUS:  IF  QR  DOESN’T  CONVERGE  IN  30 
ITERATION,  THE  PROGRAM  QUITS  AND  RETURNS  I ERR*=- 1 . 

INTERNAL  VARIABLES 

COMPLEX  X, Y, Z , El  1  , El  2 , SHIFT, ZERO 
REAL  T1 ,T2 ,T3 , AB , HMAX , HALF 
INTEGER  I , J ,K,M, ITS 

DEFINE  AB ( X )  :=  |REX|+|IMX| 

AB ( X )=ABS ( REAL ( X ) )+ABS( IMAG(X) ) 

ZERO=CMPLX ( 0.0) 

HALF=0  5 
M=IGH 
ITS=0 

0  IF(M.LE.LCW)  GOTO  300 
DO  100  K=M, LOWf  1 , -  1 

LOOK  FOR  SINGLE  SMALL  SUB- DA IGONAL  ELEMENT 

Tl=MIN( AB(H(K- 1 , K- 1 ) ) , AB ( H (K , K ) ) ) 

T3=AB(H(K,K-1 ) ) 

T2=T 1+T3 

I F ( T 1  EQ.T2 )  THEN 

HMAX=MAX(T3 , AB( H( K- 1 , K ) ) , AB ( H( K- 1 , K- 1 ) ) , AB ( H ( K , K) ) ) 

Tl=( AB(H(K- 1 ,K- 1 ) -H(K,K) ) /HMAX) • »2 
T3=M*(T3/HMAX)*(AB(H(K-1  ,K))/HMAX) 

T2=Tl+T3 

I F ( Tl  EQ  T2  )  GOTO  110 
T1=AB(H(K- l , K- 1 )+H(K,K) ) /HMAX 
T2=*T  1  +  SQRT  ( T3  ) 

I F (Tl  EQ  T2 )  GOTO  110 
END  IF 

00  CONTINUE 
10  I F ( K  EQ  M)  THEN 

ONE  EIGENVALUE  FOUND 
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GOTO  1  0 
END  1  F 

I F ( ITS  EQ. 30)  THEN 
1  ERft=-  1 
GOTO  300 
END  IF 

I F ( ITS . EQ. 1 0  OR . ITS  EQ. 20 )  THEN 
C 

C  FORM  EXCEPTIONAL  SHIFT 
C 

SHIFT=H(M,M- 1  ) 

ELSE 

C 

C  FORM  WILKINSON’S  SHIFT 
C 

Z=H(M,M-  1  )  *H(M-1  ,M) 

X=SQRT { (H(M- 1  ,M-  1 ) -H(M,M)  ) • *2+4 . 0 • Z ) 

Y=H(M-  1  , M-  I  )  +  H(M,M) 

El  1=(X+Y ) ‘HALF 

I F ( AB ( X+Y )  . LT  AB(Y-X) )  E I 1=( Y- X ) • HALF 
E I  2 _ Y  -Ell 

1F(4»AB(B12) .LT.AB(Y))  El 2=(H(M,M)  »H(M- 1 ,M- 1 )-Z)/EIl 
SHIFT=EI 1 

I F ( AB (El  1 -H(M,M)  )  . GT . A8( El  2  - H(M,M)  ) )  SHIFT=EI2 
END  IF 
ITS=ITS+1 
DO  120  1=K,M 

120  H( I , I )=H ( I , I J-SHIFT 
C 

C  QR  ITERATION  ROW  MOD  I F I  CAT I ON 
C 

DO  140  I=K,M-1 

T1=SQRT( REAL ( H( 1,1) *CONJG(H( I , I ) ) )  + 

X  REAL ( H ( 1+1,1 ) *CON JG( H ( 1+1,1)))) 

W(  I  )=H  (  I  ,  I  )  / T 1 
V( I )=H( 1+1,1) /Tl 
H( I , I )=CMPLX( Tl ) 

H(  1+1  ,  I  )=ZERO 
DO  130  J=I  +  1  ,  N 
X=H (  I  ,  J) 

H (  I  , J )=CONJG(W(  I ) ) *X+CONJG( V ( I  ) ) • H(  1+ 1  , J ) 

H(  1+1  ,  J )=W(  I)*H(l+t,J)-V(I)«X 
130  CONTINUE 

140  CONTINUE 
C 

C  COLUMN  MODIFICATION 

C 

DO  180  J=K, M-l 
DO  160  1=1  , J 
X=H ( I , J ) 

H(  l  ,  J)=*V(  J  )  •  X+V (  J  )  *H(  I  ,  J+l  ) 

H (  I  ,  J+  1  )=CONJG(W(  J) ) »H( I  ,  J+l ) -CO NJG(V( J ) ) »X 
160  CONTINUE 

H(  J+l  ,  J  )=V(  J)  *H(  J+l  ,  J+l  ) 

H( J+l  , J+l )=CONJG(W(  J ) ) *H( J+l  , J+ 1 ) 

C 

C  ACCUMULATE  TRANSFORMATION 

C 

DO  1  TO  1=1  ,N 
X=P(  I  , J  ) 

P (  I  , J )=W(  J ) *X+V( J ) • P (  I  ,  J+l  ) 

P(  I  , J+l )=CONJG(W(  J) )*P( I  , J+l  ) -  CON  J  G ( V ( J ) ) • X 
1 7  0  CONT I NUE 

180  CONTINUE 


O  O  O  M  o  o  o 
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SHIFT  BACK 

DO  210  I=K,M 
10  H( I  ,  I )=H( I  ,  I  )  +  SH I  FT 
GOTO  10 

CLEAR  THE  LOWER  TRIANGULAR  PART 

300  DO  3 SO  J=I  ,N-  1 
DO  350  t=J+l  ,N 
350  H( l , J )=ZERO 


CORDER 


1.  PURPOSE. 

The  Fortran  77  su!»outine  CORDER  re-orders  the  diagonal  of  a  complex  triangular 
matrix  S  according  to  the  real  parts  of  tau*S(i,i),  using  unitary  similarity  transforma¬ 
tions  that  preserves  triangularity.  The  transformations  are  accumulated  in  P. 


2.  USAGE. 

(A).  Calling  Sequence. 

SUBROUTINE  CORDER(S,P,tau,w,n,ns,np) 

Parameters  : 

S  is  a  complex  two-dimensional  variable  with  row  dimension  ns  and 
column  dimension  at  least  n.  On  input,  S  contains  a  complex  triangular 
matrix  of  order  n  (a  Schur  form  of  B).  On  output,  S  contains  the  tri¬ 
angular  matrix  which  is  unitarily  similar  to  the  previous  one  and  has  its 
diagonal  elements  ordered  according  the  the  real  parts  of  the  diagonal 
of  tau*S,  where  tau  is  a  complex  parameter. 

P  is  a  complex  two-dimensional  variable  with  row  dimension  np  and 
column  dimension  at  least  n.  On  input,  P  contains  an  unitary  matrix 
that  transform  B  to  its  Schur  form  S:  S=»P*BP.  On  output,  P  con¬ 
tains  the  updated  unitary  matrix  that  transform  B  to  the  new  S. 

w  is  a  real  one-dimensional  variable  of  order  n.  Work  space. 

tau  is  an  input  complex  number,  parameter. 

n  is  an  input  integer  equal  to  the  order  of  the  matrix  B. 

ns,np  are  input  integers  equal  to  the  row  dimension  of  S  and  P  respectively. 


S.  DISCUSSION  OF  METHOD  AND  ALGORITHM. 


The  re-ordering  of  the  diagonals  of  S  accomplished  by  a  sequence  of  adjacent  swaps.  A 
bubble  sort  is  performed  on  elements  w(i):=Re(tau*S(i,i)).  Whenever  a  swap  is 
needed,  say,  to  swap  w(k)  and  w(k+  1),  an  unitary  matrix  U  is  constructed  so  that 


U*- 


sll 

s  1 2 

s22 

s  1 2 

0 

s22  U  = 

0 

sll 

here  s  1 1 ,  sl2  and  s22  are  the  current  elements  S(k,k),  S(k,k+  1)  and  S(k+  l,k+  1) 
respectively.  It  is  elementary  to  show  that  the  U  defined  below  effects  such  a  transfor¬ 
mation: 
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U 


if  sll=s22  ; 


if  sl2=0  and  sll^s22  ; 


x  -x/x 

1  X 


x  := 


s!2 

s22-sll 


otherwise. 


To  complete  the  transformation,  columns  k  and  (k+  1)  of  both  S  and  P  are  postmulti- 
plied  by  U;  rows  k  and  (k+  1)  of  the  new  S  are  premultiplied  by  U*. 

The  algorithm  terminates  when  the  bubble  sort  is  finished,  or  when  no  swap  occurs  in 
any  sweep  (which  means  that  all  the  elements  are  in  order). 

There  are  approximately  67  executable  Fortran  statements. 


4.  REFERENCES. 
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5.  PROGRAM  LISTING. 


SUBROUT INE  CORDER ( S , P , TAU , W,  N , NS ,NP) 

COMPLEX  S(NS ,N) . P(NP,N) , TAU 
REAL  \V(  N ) 

I NTEGEK  N , NS , NP 

THIS  SUBROUTINE  RE-ORDERS  THE  DIAGONAL  OF  A  COMPLEX  TRIANGULAR  MATRIX  S 
BY  UNITARY  TRANSFORMATIONS  (KEEPING  S  TRIANGULAR).  THE  NEW  ORDERING  IS 
ACCORDING  TO  THE  REAL  PARTS  OF  TAU«T( 1,1)  THE  TRANSFORMATIONS  ARE 
ACCUMULATED  IN  P 

CODE  IN  F7  7  BY  K.C  NG ,  UC  BERKELEY,  REVISED  ON  5/20/85 

F7  7  GENERIC  FUNCTIONS  ABS  ,  CONJG ,  REAL,  SQRT 

F77  NON-GENERIC  FUNCTION;  CMPLX 

GLOSSARY: 

S  AN  UPPER  TRIANGULAR  COMPLEX  MATRIX  WHOSE  DIAGONALS  ARE  GOING 

TO  BE  RE-ORDERED 

P  A  UNITARY  MATRIX  IN  WHICH  THE  TRANSFORMATIONS  ARE  ACCUMULATED 

TAU  COMPLEX  PARAMETER 

W  COMPLEX  ARRAYS.  WORK  SPACE 

N  DIMENSION  OF  MATRIX  S.P 

NS , NP  LEADING  DIMENSION  OF  ARRAYS  S  AND  P  RESPECTIVELY 

INTERNAL  VARIABLES 

I NTEGER  I , J , K , NSWAP 
COMPLEX  X.Y.U21 ,Ull ,UI2, ZERO 
REAL  RX , RY , SQ , ONE 
2ERO=CMPLX( 00) 

ONE=I 

DO  10  1=1  ,N 

0  W(  I  )=REAL ( TAU • S ( I  ,  I  ) ) 

DO  1 10  1=2 , N 

NSWAP=0 

DO  10  0  It=N  ,1,-1 
IF  (W(K)  .  LT  W(K-  1  )  )  THEN 
NSWAP=NSWAP+  1 
RX=W(K) 

W(K)=W(K-  1  ) 

W(K- 1 )=RX 

•••••  CONSTRUCT  THE  UNITARY  MATRIX  U  ••••• 

X=S(K- 1 , K ) 

Y=S(K,K) - S ( K- 1 , K- 1 ) 

IF  (X  EQ  ZERO)  THEN 
Ul  1=0  0 
U2  1=1  0 
Ul  2=U2  1 
ELSE 

RX=ABS ( X ) 

RY=ABS(Y) 

IF  (RX  LE  RY)  THEN 

SQ=l  0+ ( RX/RY ) • • 2 

IF  (  SQ  NE  ONE )  GOTO  3 0 

X  <  1 Y  AND  I X  Y i  «  1  (  I  E  1  *  X  Y ,  2  =  1  ) 

Ul  1  =X  /  Y 
U2  1=1  0 

Ul 2=- ( X, CONJG(X) ) • ( CONJG ( Y ) , Y ) 


c  o  r  d  e  r 


SQ=1  .  0+  ( RY/RX  )  •  *  2 
IF  (SQ.NE.ONE)  GOTO  30 

I X  I  >  | V  I  AND  | X/Y |  »  1  (I  E..  1+ (X/ Y | *  2  =  |X/Y|'2) 

U1 1=(X/CMPLX(RX) ) • (CMPLX(RY)/Y) 

U2 1=CMPLX( RY/RX) 

U 1 2=-  U 1  1  •  CON  J  G ( Y / X ) 

END  IF 
END  IF 
GOTO  40 

1 X/ Y |  IS  NOT  TOO  BIG  AND  NOT  TOO  SMALL 

X=X/Y 

U2 1=CMPLX( 1.0/ SQRT ( 1  0+ ( RX/RY ) • • 2 ) ) 

Ul 1=X*  U2 1 

Ul  2=- U 1 1 /CON JG ( X ) 


END  OF  THE  CONSTRUCTION  OF  U 


C  COLUMN  MODIFICATION 
C 

4  0  DO  6  0  J= 1  . K- 2 

X  =S  ( J  ,K-  1 ) *U 1 1+S ( J , K) *U2 1 

S  (  J  ,  K )  ==S  (  J  ,K-  1  )  •  U 1 2+S (  J , K) • U 1 1 

S ( J , K- 1 )=X 
60  CONTINUE 

C 

C  ACCUMULATE  THE  TRANSFORMATION 
C 

DO  7  0  J=1  , N 

X  =P( J ,K- 1 ) *UI 1+P( J ,K) *U21 

P ( J , K )  =P( J.K-1 )  *Ul  2+P( J , K) *Ul 1 

P(  J ,K-  1  )=X 
70  CONTINUE 

C 

C  RCW  MOD  I F I  CAT  I  ON 

C 

Ul l=CON  JG ( U 1 1  ) 

U2 l=CONJG( U2 1 ) 

U 1 2=CON  JG ( U 1 2 ) 

DO  8  0  J=K+ 1  . N 

X  =S(K- 1 , J ) *UI 1+S(K. J ) *U2 1 

S  (  K  ,  J  )  =S  (  K  -  1  .  J  )  •  U  1  2+  S  (  K  .  J  )  •  U  1  1 
S  (  K-  1  .  J  )=X 
80  CONTINUE 

C 

C  SWAP  THE  DIAGONAL 
C 

X  =S  ( K  .  K ) 

S ( K , K )  =S(K- 1 ,K- l ) 

S  (  K-  1  ,  K-  1  )=X 
END  IF 

100  CONTINUE 

IF  ( NSWAP  EQ  0)  RETURN 
110  CONTINUE 
RETURN 
END 
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CEXPHY 


1.  PURPOSE. 

The  Fortran  77  subroutine  CEXPHY  computes  the  exponential  of  a  complex  triangular 
matrix  E=exp(tau*B)  by  first  determining  a  block  diagonal  clustering  on  B,  calling 
subroutine  CEXPRI  to  compute  the  exponential  of  the  diagonal  blocks  and  finally 
applying  Parlett’s  recurrence  formula  to  fill  up  the  all  the  off-diagonal  elements. 


2.  USAGE. 

Calling  Sequence. 

SUBROUTINE  CEXPHY(B,E,w,ip,idx,tau,ovft,ierr,n,nb,ne) 

Parameters  : 

B  is  an  input  complex  two-dimensional  variable  with  row  dimension  nb 
and  column  dimension  at  least  n.  B  contains  a  complex  triangular 
matrix  of  order  n,  with  its  diagonal  elements  ordered  according  to  the 
real  part  of  tau*B(i,i)  (cf.  subroutine  CORDER). 

E  is  a  complex  two-dimensional  variable  with  row  dimension  ne  and 
column  dimension  at  least  n.  On  output,  E  contains  the  exponential  of 
tau*B. 

w  is  a  complex  one-dimensional  variable  with  dimension  at  least  on. 

*  Work  space. 

ip.idx  are  integer  arrays  of  dimension  n.  Work  space. 

ovft  is  a  real  number  equal  to  the  overflow  threshold. 

tau  is  an  input  complex  number,  parameter. 

ierr  is  an  integer  variable  indicating  error  condition.  If  the  exponential  of 
some  eigenvalues  overflow,  ierr  is  set  equal  to  -2,  and  the  subroutine 
terminates. 

n  is  an  input  integer  equal  to  the  order  of  the  matrix  B. 

nb,ne  are  input  integers  equal  to  the  row  dimension  of  B  and  E  respectively. 
Subprograms  :  CEXPRI 

S.  DISCUSSION  OF  METHOD  AND  ALGORITHM. 

Let  z(l),z(2) . z(n)  denote  the  eigenvalues  of  tau*B.  CEXPHY  first  determines  a  clus¬ 

tering  of  the  diagonals  of  E=exp(tau*B)  (see  Fig.  I  in  the  next  page).  The  requirement 
is  that  z(i)  in  different  clusters  must  be  well  separated.  A  typical  criterion  reads  as  fol¬ 
lows  :  the  i-th  and  j-th  diagonal  elements  of  E  are  in  the  same  cluster  if 


|z(iH(j)l  <  g(H)  • 

We  have  found  that  a  second  degree  polynomial  is  sufficient  :  g(k)=aj+  a^k2. 
(See  the  program  listing  for  the  values  of  the  constants  at.)  This  criterion  arose  from 
studying  the  computation  of  the  exponential  divided  differences.  (In  [2]  it  suggests  g(k) 
g(k)=k*(l+  0.1*ln  k);  however,  we  found  that  a  quadratic  polynomial  is  more  realis¬ 
tic.) 

If  two  clusters  so  determined  overlap,  then  we  merge  them  together.  The  following 
figure  is  a  typical  clustering  of  E. 


XX 

X . 

XXX 
XX  - 
X  ■  ■ 

X  X 
X 


X  computed  by  CEXPNT 

computed  by  Parlett's  recurrence 


Fig.  1.  Computing  exp(tau*B)  by  hybrid  method. 


After  the  clusters  are  determined,  the  subroutine  CEXPRI  is  called  to  compute  the 
diagonal  blocks  of  E  (see  the  description  of  CEXPRI  for  details  of  the  method). 

Finally,  the  rest  of  E  is  computed  by  Parlett's  recurrence  (cf.  [3]  or  [l])  : 

B(i,jHE(j,j)-E(i,i))+  £  (B(i,k)*E(k,j)  -  E(i,k)*B(k,j)) 

'  (B(j,j)-B(i,i)) 

The  recurrence  is  obtained  from  the  commutativity  of  E  and  B.  (For  real  quasi  tri¬ 
angular  B,  the  diagonal  may  contain  some  2x2  blocks.  In  that  case,  we  solve  BE- 
EB=0  directly  for  each  block  in  turn.) 

There  are  approximately  55  executable  Fortran  statements. 


I  REFERENCES. 

(1}  G.H.  Golub  &  C.F.  VanLoan,  Matrix  Computations,  Johns  Hopkins  University 
Press,  1983. 

(2)  A.  McCurdy,  K.C.  Ng  and  B.N.  Parlett,  Accurate  Computation  of  Divided 
Differences  of  The  Exponential  Function,  Mathematics  of  Computation, 
Volume  43,  Number  168,  October  1984,  pp501-528. 


[3]  B.N.  Parlett,  A  Recurrence  Among  the  Elements  of  Functions  of  Triangular 
Matrices,  Linear  Algebra  Appl.,  14(1976),  pp. 117-121. 
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5.  PROGRAM  LISTING. 


SUBROUTINE  CEXPHY ( B , E ,  W,  IP,  I DX , TAU , OVFT ,  I  ERR , N , NB , NE ) 

COMPLEX  B(NB,N),  E(NE,N),  W(S*N),  TAU 
REAL  OVFT 

INTEGER  N.NB.NE, i ERR , IDX(N) . IP(N) 

THIS  SUBROUTINE  COMPUTES  THE  EXPONENTIAL  OF  (COMPLEX  TRIANGULAR)  B  BY 
CEXPRI  AND  PARLETT'S  RECURRENCE  RESULT  STORED  IN  E 
CODE  IN  F77  BY  K  C  NG ,  UC  BERKELEY,  REVISED  ON  5/ 20/85 

F7  7  GENERIC  FUNCTIONS  ABS  ,  I MAG.  REAL,  MAX 

REQUIRED  SUBROUTINE  CEXPRI 

LOCAL  REAL  VALUED  FUNCTIONS  AB  G 

GLOSSARY : 

B  AN  UPPER  TRIANGULAR  COMPLEX  MATRIX 

E  ON  OUTPUT,  THE  EXPONENTIAL  OF  B 

W  COMPLEX  ARRAYS ,  WORK  SPACE 

IP , IDX INTEGER  ARRAYS,  WORK  SPACE 
TAU  COMPLEX  PARAMETER 

OVFT  COMPUTER  (REAL)  OVERFLOW  THRESHOLD 

I  ERR  (INTEGER)  ERROR  STATUS  RETURN  I ERR=  -2  IF  SOME  EIGENVALUE 
IS  TOO  LARGE  TO  EXPOENET I  ATE 
N  DIMENSION  OF  MATRIX  E 

NB , NE  LEADING  DIMENSION  OF  ARRAYS  B  AND  E  RESPECTIVELY 
INTERNAL  VARIABLES 
COMPLEX  Z 

INTEGER  I . J  .  I  1  .  J  1  ,  K 
REAL  G.A1  ,*A2  .A3  ,  AB  ,  ZERO 

DATA  A 1  , A2 , A3 , ZERO/ -1797804884,1.958414640,  02930024390,0  0/ 

DEFINE  REAL  VALUED  FUNCTION  G  AND  AB 

G(  I  )=(  A1+  I  •  (  A2+  I  •  A3  )  ) 

AB  (  Z  )=ABS  (  REAL  (  Z  )  )  +  ABS(  IMAG(Z)  ) 

INITIALIZE  (SAVE  THE  EIGENVALUES  IN  THE  LAST  COLUMN  OF  E  TEMPORARILY) 

DO  2  0  J=1  ,  N-  1 
DO  20  1=1,1 

0  E( I , J )=6( I , J ) *TAU 
DO  4  0  1=1  ,N 
0  E ( I , N)=6 ( 1,1) »TAU 

MAKE  ROOM  IN  B 

DO  100  J=1  ,N/2 
K=N-  J 

DO  8  0  1=1  ,  J 

0  B(K-t-I  ,K)=B(  I  ,  J) 

00  CONTINUE 

FIND  THE  NEXT  CLUSTER  (Il.Jl) 


120 


C 


1  1=1 
J  1  =1  t 

I F (  I  1  . GT  N )  GOTO  2  6  0 
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c  ...  FOR  1=11,  1+1  WHILE  I  <  3  1+1  DO  •••• 

C 

1  =  1  1 

MO  I F  (  1  .  LE  J  1  )  THEN 

DO  16  0  3=N , 3 1  , -  1 

IF(AB(£( I ,N) -E( J.N) )  LE  G(3-I))  GOTO  180 
160  CONTINUE 

180  3l=MAX(  3 , 31 ) 

IF( 3 1 . EQ  N)  GOTO  200 
1=1+1 
GOTO  1 4  0 
END  IF 
C 

C  COMPUTE  THE  EXP  OF  THE  CLUSTER 
C 

200  CONTINUE 

I F ( 3 1 . EQ . N )  THEN 
DO  220  1  =  1  , N 

220  E( I ,N)=B( I , N ) *TAU 

END  IF 

CALL  CEXPRI ( E( I  1 .  I  1  )  ,W, I P ( I  1 ) ,  I DX( I  1 )  , B.OVFT , I  ERR , J 1 
DO  240  1=1 1 , J 1 

240  I DX ( I )=l 1 
I  1=J  1+1 
GOTO  120 
260  CONTINUE 
C 

C  RESTORE  B 
C 

DO  300  J=1 ,N/2 
K=N-  J 

DO  2  80  1=1  ,  3 

280  B ( I , J )=B(K+ I ,K) 

300  CONTINUE 

DO  3  20  1=1  ,N 
320  W(  I  )=TAU*B(  1,1) 

C 

C  FILL  UP  THE  OFF-DIAGONAL  BLOCK  BY  RECURRENCE 
C 

DO  380  I=N, 1,-1 
DO  380  3=1+1 ,N  . 

I F ( I DX ( I  ) - NE .  1DX( J )  OR  AB(W(  J  )  -W(  I ) )  GT  G( 3 -  I ) ) 
Z=B( I , J) *(E( J , J )-E( I , I ) ) 

DO  340  K=I+ 1 .3-1 

Z=Z+B( I ,K)»E(K,J)-E( I ,K)*B(K,J) 

340  CONTINUE 

Z=Z/(B( J, J } - B ( 1,1)) 

E( I  , J  )=Z 
END  IF 

360  CONT I NUE 
380  CONTINUE 
C 

C  CLEAR  THE  LOWER  B  AND  E 
C 

DO  440  3=1  ,N-  1 

DO  4  00  I=J+1  , N 

400  E( I , J )=ZERO 

DO  4  20  I=J+1  , N 
420  B( I , J )=ZERO 

440  CONTINUE 
RETURN 
END 


-  I  1+1  ,NE) 


THEN 
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CEXPRI 

(Subprogram  of  CEXPHY) 


1.  PURPOSE. 

The  Fortran  77  subroutine  CEXPRI  computes  the  exponential  of  a  complex  triangular 
matrix  E  using  a  combination  of  the  Newton  polynomial  method  and  Parlett’s 
recurrence  formula.  The  result  overwrites  E. 


2.  USAGE. 

Calling  Sequence. 

SUBROUTINE  CEXPRI(E,w,ip,idx,s,ovft,ierr,n,ne) 


Parameters  : 

E  is  a  complex  two-dimensional  variable  with  row  dimension  nb  and 
column  dimension  at  least  n.  On  input,  E  contains  a  complex  triangular 
matrix  of  order  n.  On  output,  exp(E)  overwrites  E. 

w  is  a  complex  one-dimensional  variable  with  dimension  at  least  5n. 
Work  space. 

ip  is  an  integer  array  of  dimension  n.  Work  space, 

idx  is  an  integer  array  of  dimension  n.^Work  space. 

s  is  an  one  dimensional  complex  array  of  dimension  at  least  n(n-2)/4. 

Work  space. 

ovft  is  a  real  variable  set  equal  to  the  overflow  threshold. 

ierr  is  an  integer  variable  indicating  error  condition.  If  the  exponential  of 
some  eigenvalues  overflow,  ierr  is  set  equal  to  -2,  and  the  subroutine 
terminates. 

n  is  an  input  integer  set  equal  to  the  column  dimension  of  the  matrix  E. 

ne  is  an  input  integer  set  equal  to  the  row  dimension  of  E  respectively. 

Subprograms  :  CEXPNT,  CINDEX,  CMSWAP 

S.  DISCUSSION  OF  METHOD  AND  ALGORITHM. 

Let  z(l),z(2),...,z(n)  denote  the  eigenvalues  of  E.  CEXPRI  first  determines  a  grouping 
of  z(i)  by  subroutine  CINDEX.  The  criterion  is  :  if  the  imaginary  parts  of  z(i)  and  z(j) 
are  differ  by  less  than  pi*3. 14159...  then  z(i)  and  z(j)  should  belong  to  the  same  group. 
We  use  an  integer  array  idx(i)  to  index  z(i).  Thus,  if  z(i)  and  z(j)  should  be  together, 
we  set  idx(i)=»idx(j). 
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After  determining  idx(i),  we  apply  CMSWAP  if  necessary  to  have  E’s  diagonal  ordered 
according  to  idx(i)  (via  untary  similarity  transformations).  Then  on  each  group  (diago¬ 
nal  block)  subroutine  CEXPNT  is  called  to  compute  its  exponential.  Off  diagonal  ele¬ 
ments  are  filled  up  by  Parlett’s  recurrence  (cf.  [3])  (resutls  overwrite  E). 

Finally,  we  undo  the  unitary  similarity  transformations  on  E. 

There  are  approximately  42  executable  Fortran  statements  in  CEXPRI,  31  in  subrou¬ 
tine  CMSWAP,  and  42  in  subroutine  CINDEX. 


4.  REFERENCES. 

[1]  G.H.  Golub  &  C.F.  VanLoan,  Matrix  Computations,  Johns  Hopkins  University 
Press,  1983. 

(2]  A.  McCurdy,  K.C.  Ng  and  B.N.  Parlett,  Accurate  Computation  of  Divided 

Differences  of  The  Exponential  Function,  Mathematics  of  Computation, 
Volume  43,  Number  168,  October  1984,  pp501-528. 


[3]  B.N.  Parlett,  A  Recurrence  Among  the  Elements  of  Functions  of  Triangular 
Matrices,  Linear  Algebra  Appl.,  14(1976),  pp. 117-121. 
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5.  PROGRAM  LISTING. 


SUBROUTINE  CEXPRI (E,W,  IP,  IDX.S.OVFT, IERR.N.NE) 

COMPLEX  E(NE,N)  , S (N» (N- 2 ) / 4 )  ,W(S*N) 

INTEGER  IP(N) ,  IDX(N)  ,N,NE,  I  ERR 
REAL  OVFT 

THIS  SUBROUTINE  COMPUTES  THE  EXPONENTIAL  OF  (COMPLEX  TRIANGULAR)  E  BY 
CEXPNT  AND  PARLETT  RECURRENCE  RESULT  OVERWRITES  E. 

CODE  IN  F77  BY  K.C  NG ,  UC  BERKELEY;  REVISED  ON  6/6/85. 

F77  GENERIC  FUNCTIONS  CONJG 

REQUIRED  SUBROUTINES  CINDEX,  CEXPNT,  CMSWAP 

GLOSSARY : 

E  ON  INPUT,  AN  UPPER  TRIANGULAR  (COMPLEX)  MATRIX,  ON  OUTPUT, 

THE  EXPONENTIAL  OF  E 
W,S  COMPLEX  ARRAYS,  WORK  SPACE 
IP, I DX INTEGER  ARRAYS,  WORK  SPACE 
OVFT  COMPUTER  (REAL)  OVERFLOW  THRESHOLD 

I  ERR  (INTEGER)  ERROR  STATUS:  RETURN  I ERft=  -2  IF  SOME  EIGENVALUE 
IS  TOO  LARGE  TO  EXPONENTIATE. 

N  DIMENSION  OF  MATRIX  E 

NE  LEADING  DIMENSION  OF  ARRAY  E 

INTERNAL  VARIABLES: 

INTEGER  I.J.K.L 
COMPLEX  X 
DO  2  0  1=1  ,  N 

0  W(  I  )==£(  1,1) 

FIND  OUT  THE  NEW  ORDERING  FOR  W(  I )  » 

CALL  CINDEX(N,W,W(NM  )  , W(  2 »N+ 1 ) , I P ,  I DX ) 

SWAP  E( I , I )  ACCORDING  TO  THE  NEW  ORDERING 


L=0 

DO  100  I=N, 1,-1 

DO  4  0  J=I  ,1,-1 

l F ( I DX ( J ) . EQ  I P ( I ) )  GOTO  6  0 
0  CONTINUE 

0  DO  8  0  K=J+1  ,  I 

X=E(K- 1 ,K) / ( E(K,K) -E(K- 1 ,K- 1 ) ) 

L=L+  1 

S(L )=CONJG(X) 

CALL  CMSWAP (E,W,  I DX , X , K ,N , NE ) 

0  CONTINUE 

1P( I )=J+1 
00  CONTINUE 

DO  110  1=1 ,N- 1 
DO  110  J=I+  1  ,  N 
10  E( J , I )=£( 1  .  J) 

COMPUTE  THE  EXPONENTIAL  OF  THE  DIAGONAL  BLOCKS 

K=IDX( 1 ) 

1  =  1 

DO  140  J=1  ,  N 

TF( J . EQ  N)  THEN 

CALL  CEXPNT ( E( 1,1) ,W(  I )  ,OVFT,  IERR.N- 1+1  , NE ) 
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ELSE  I F ( 1 DX ( J )  NE  K)  THEN 

CALL  CEXPNT ( E ( 1  , I ) ,W(  I ) ,OVFT,  I  ERR , J- 1  , NE) 
l=J 

K=IDX(  I ) 

END  IF 

I F ( i ERR . EQ  -2 )  RETURN 
140  CONTINUE 
C 

C  COMPUTE  THE  OFF-DIAGONAL  ELEMENTS  BY  PARLETT  RECURRENCE 
C 

DO  300  I=N,  1,-1 
DO  300  J— 1+  1  , N 

I F ( I DX ( I )  NE  IDX(J))  THEN 
X=£  ( J ,  I ) • ( E ( J , J ) - E (  1  .  I ) ) 

DO  2  80  10=1+1  ,  J-  1 

280  X=X+E(K, I ) «  E (K , J) -E( I ,K) •  E( J ,K) 

X=X/  ( W(  J  )  -W(  1  )  ) 

E( 1 , J )=X 
END  I  F 

300  CONTINUE 
C 

C  SWAP  BACK 
C 

DO  320  1=1 ,N 
DO  3  2  0  K=I  ,  I P  (  I  ) ,  -  1 
'  X=S  (  L  ) 

L=L  -  1 

CALL  CMSWAP(E,W,  IDX.X,  -K.N.NE) 

320  CONTINUE 
RETURN 
END 


‘SUBROUTINE  CMSWAP(E,Z  ,  I  DX  ,  X  ,  K ,  N ,  NE ) 

COMPLEX  E(NE.N) , Z(N) ,X 
INTEGER  K.NE.N, IDX(N) 

C 

C  THIS  SUBROUTINE  SWAPS  THE  K  AND  K- 1  DIAGONAL  ELEMENTS  OF  MATRIX  E  USING 
C  UNITARY  C  TRANSFORMATION  (THE  TRANSFORMATION  IS  ENCODED  IN  THE  INPUT  X) 
C  IT  ALSO  SWAPS  THE  CORRESPONDING  ELEMENTS  IN  ARRAY  Z  AND  ARRAY  IDX, 

C  CODE  IN  F77  BY  K  C  NG ,  UC  BERKELEY;  REVISED  ON  5/18/85. 

C 

C  F77  GENERIC  FUNCTIONS  ABS ,  CONJG,  REAL,  SQRT 

C  F77  NON- GENERIC  FUNCTION:  CMPLX 

C 

C  GLOSSARY 

C  E  AN  UPPER  TRIANGULAR  (COMPLEX)  MATRIX 

C  Z  COMPLEX  ARRAY 

C  IDX  INTEGER  ARRAY 

C  X  COMPLEX  NUMBER  CONTAINS  THE  TRANSFORMATION  INFORMATION 

C  K  INTEGER  INDEX 

C  N  DIMENSION  OF  MATRIX  E 

C  NE  LEADING  DIMENSION  OF  ARRAY  E 

C 

C  INTERNAL  VARIABLES 
C 

COMPLEX  U1I .U12.U21 .ZERO 

I  NT EGER  I  , J 

DATA  ZERO / (00,00)/ 

C 

C  RECONSTRUCT  THE  UNITARY  MATRIX  FROM  X 
C 

U2 1=CMPLX( l  0 / SQRT ( 1  0+REAL ( CONJG (X ) *X) ) ) 


o  o  o  o  o  o  o  o  o  o  o  o  o  o  a  o  o  ooo«  ooo 


c  exphy - c  exp  r  i 


Ul 1=X*U21 
Ul2=U2  1 

1 F ( X . NE . ZERO)  THEN 
IF(K.GT.O)  THEN 
U 1 2=- U 1 1 /CONJG(X) 

ELSE 

U2  1=- U 1 1 / CONJG( X ) 

END  IF 
END  IF 
K=*ABS(K) 

PERFORM  E«U 

DO  3  0  1=1  ,  K-  2 

X  =E(  I  ,K- 1 ) *U1  1+E(  1 ,K) *U2 1 

E( I , K)  =E( I ,K- 1 ) *U12+E( I ,K) *U1 1 
E(  I , K- 1 )=X 
0  CONTINUE 

PERFORM  UH*  E 

Ul l=CONJG( U 1  1  ) 

U2  l=CONJG(  U2I  ) 

U 1 2=CONJG( Ul 2 ) 

DO  4  0  J=K+ 1  , N 

X  =E(K- 1 , J ) • Ul 1+E(K, J )*U21 

E(K,  J )  =E(K- 1 ,  J ) »U1 2+E(K,  J ) • U 1 1 

E(K- 1 , J  )=X 
0  CONTINUE 

SWAP  THE  DIAGONAL  OF  E  AND  THE  CORRESPONDING  ELEMENTS  IN  Z  AND  1 DX 


I 

=IDX(K) 

IDX(K) 

=IDX(K- 1 ) 

I DX(K- 1 ) 

=1 

X 

=E(K,K) 

E(K,K) 

=E(K- 1 ,K- 1 

E(K- 1 , K- 1 

)=x 

X 

=Z(K) 

Z  (K) 

=Z  { K-  1  ) 

Z(K-l) 

=X 

RETURN 

END 

SUBROUTINE  CINDEX(N,W,  Z I  ,  IQ,  IP.  IDX) 

COMPLEX  W(N) 

REAL  Z I ( N) 

INTEGER  IPiN) , IQ(N) , IDX(N) ,N 

THIS  SUBROUTINE  ASSIGNS  A  GROUP  NUMBER  IDX(I)  TO  EACH  OF  THE  IMAGINARY 
PARTS  Z I  (  I  )  OF  THE  EIGENVALUES  W(  I )  THE  CRITERION  IS  IF  THE  DISTANT 
BETWEEN  ZI(I)  AND  Z I ( J )  ARE  LESS  THEN  P!=3  14158  THEN  THEY  HAVE 
THE  SAME  GROUP  NUMBER 

CODE  IN  F77  BY  K  C  NG ,  UC  BERKELEY,  REVISED  ON  5/18/85 
F7  7  GENERIC  FUNCTIONS  ABS  ,  ACOS ,  IMAG 

GLOSSARY 

W  INPUT  COMPLEX  ARRAY  CONTAINING  THE  EIGENVALUES 

ZI  REAL  ARRAY  CONTAINING  THE  IMAGINARY  PARTS  OF  W 
IDX  OUTPUT  INTEGER  ARRAY  CONTAINING  THE  GROUP  NUMBERS  OF  THE 
CORRESPONDING  ELEMENTS  IN  ZI 
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C  IP  OUTPUT  INTEGER  ARRAY  CONTAINING  THE  NEW  ORDERING  OF  ZI 

C  IQ  INTEGER  ARRAY,  WORK  SPACE 

C  N  DIMENSION  OF  ARRAY  W,  Z I  , I P , I Q 

C 

C.  INTERNAL  VARIABLES 
C 

REAL  PI , ZERO 

INTEGER  I , J ,K,L , IDIFF, ITEMP 
DATA  ZERO/0 . 0/ 

P 1=2  . 0*ACOS(ZERO) 

DO  10  1=1 , N 

I Q ( I )=■+! 

ZI  (  I  )  =  1MAG(W(  I  )  ) 

10  CONTINUE 

!P( 1  )=* 

I DX( 1 )=1 
L=1 
K=1 
1=1 
C 

C  DETERMINE  I DX ( I ) 

C 

20  IFjL.LT  N)  THEN 
30  I F ( I  LE  N-L)  THEN 

I F ( ABS ( 2  I (  I Q ( I ) ) - Z I ( I P( K) ) )  GE  PI  )  THEN 
1  =  1+1 
GOTO  30 

ELSE 

ITEMP=IDX( I P(K) ) 

END  IF 
ELSE 
1=1 
K=K+  1 

IF(K.LE.L)  GOTO  30  - 
ITEMP=!Q(I)  * 

END  IF 
L=L+1 

IP(L)=IQ( I ) 

I  DX (  I  P  (  L  )  )  =  I  TEMP 
DO  4  0  J=I  ,  N-  L 

40  iq( j )=iq( j+i ) 

GOTO  20 
END  IF 
C 

C  PUT  THE  NEW  ORDERING  IN  I P { ) 

C 

I DI FF=0 
DO  12  0  J=2  ,  N 
DO  120  1=1  , J -  1 

!F(  IDX( J  )  GT  I DX(  I ) )  THEN 
I D 1 FF=I  D I FF -  1 

ELSE  I F (  I DX ( J )  L T .  I DX (  I ) )  THEN 
I D I F  F= I D I FF+ 1 
END  IF 

120  CONTINUE 
K=1 

DO  160  L=1  , N 
J=i 

I F ( IDIFF  GT  0 )  J=N+ 1 -L 
DO  140  1=1 ,N 

I F ( I DX ( I )  EQ  J )  THEN 
I  P  (  K )  =  J 
K=K^  1 
END  IF 
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CONTINUE 

CONTINUE 

RETURN 

END 
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CEXPNT 

(Subprogram  of  CEXPRI) 

1.  PURPOSE. 

The  Fortran  77  subroutine  CEXPNT  computes  the  exponential  E  of  a  complex  tri¬ 
angular  matrix  using  the  Newton  interpolating  polynomial. 

2.  USAGE. 

(A).  Calling  Sequence. 

SUBROUTINE  CEXPNT(B,w,ovft,ierr,n,nb) 

Parameters  : 

B  is  a  complex  two-dimensional  variable  with  row  dimension  nb  and 
column  dimension  at  least  n.  On  input,  B  contains  a  complex  triangular 
matrix  of  order  n,  with  its  diagonal  elements  ordered  according  to  the 
real  part  of  B(i,i).  On  output,  it  contains  E  . 

w  is  a  complex  one-dimensional  variable  with  dimension  at  least  5n. 
Work  space. 

ovft  is  the  overflow  threshold. 

ierr  is  an  integer  variable  indicating  the  error  condition.  If  the  exponential 
of  some  eigenvalues  overflow,  ierr  is  set  equal  to  -2,  and  the  subroutine 
terminates. 

n  is  an  input  integer  equal  to  the  order  of  the  matrix  B. 

nb  is  an  input  integer  equal  to  the  row  dimension  of  B. 

Subprogram  :  CDDEXP. 

S.  DISCUSSION  OF  METHOD  AND  ALGORITHM. 

Before  computing  E=exp(B),  we  scale  down  the  size  of  B  so  that  the  the  maximum 
spread  of  the  imaginary  parts  of  B’s  eigenvalues  is  less  then  5.  Then  exp(B)  is  com¬ 
puted  by  the  Newton  polynomial  method  from  the  bottom  row  to  the  top  row. 

Suppose  that  row  j  of  exp(B)  has  been  computed.  To  compute  row  j-1,  we  first  evalu¬ 
ate  the  scaled  divided  differences  k!A,*exp  for  k=0,...,n-i  where  A,*exp  denotes  the  k-th 
divided  difference  of  exp  at  B’s  eigenvalues  z(i),z(i+  l),...,z(i+  k)  (z(i)  is  the  i-th  diago¬ 
nal  element  of  B).  T^is  is  done  by  calling  subroutine  CDDEXP.  Then  E(i,i),...,E(i,m) 
are  computed  using  the  Newton  interpolating  polynomial  as  follows. 

Let  B(  denote  the  principal  submatrix  of  B  in  rows  i  through  n,  of  B 
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z(i)  ^1,1+1  biiD 

0  z(i+ 1)  ■  •  bi+liB 

B,  =  0  0  • 

0  0  0- 

0  0  0  0  z(n) 

The  i-th  row  of  exp(B)  can  be  regarded  as  the  first  row  of  exp(Bi);  which  can  be 
represented  by 

exp(B,)  =  exp(z(i))I  +  (k!Aikexp)(-^-  JJ  (B(  -  z(j)-I)).  (1) 

lc=*l  K-  j=i 

Only  one  complex  vector  is  needed  to  accumulate  the  top  row  of  the  product  matrix 
J|(Brz(j)  I)  and  hence  the  top  row  of  exp(Bj).  The  above  is  repeated  until  i  =  1. 

Finally  any  needed  squarings  are  performed. 

There  are  approximately  115  executable  Fortran  statements. 


4.  REFERENCES. 
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5.  PROGRAM  LISTING. 


SUBROUTINE  CEXPNT ( B ,  W,  OVFT ,  I  ERR , N , NB ) 

COMPLEX  B ( NB , N ) ,  W(S*N) 

INTEGER  N, NB ,  I  ERR 
REAL  OVFT 

THIS  SUBROUTINE  USES  THE  NEWTON  POLYNOMIAL  METHOD  WITH  SCALING  AND 
SQUARING  TO  COMPUTE  THE  EXPONENTIAL  OF  B  RESULTS  OVERWRITE  B  THE 
OLD  MATRIX  B  IS  KEPT  IN  THE  LCWER  TRIANGULAR  PART  OF  B 
CODE  IN  F  7  7  BY  K  C  NG ,  UC  BERKELEY;  REVISED  ON  6/1 7 / 8  S 

F77  GENERIC  FUNCTIONS  ABS  ,SIN,  EXP,  IMAG ,  MAX,  MIN.  REAL 

F77  NON-GENERIC  FUNCTION  CMPLX 

REQUIRED  SUBROUTINES  CDDEXP 

REQUIRED  BLAS  SUBROUTINES  CAXPY ,  CCOPY ,  CSCAL ,  CSSCAL 
REQUIRED  BLAS  FUNCTION  CDOTU 

GLOSSARY : 

B  ON  INPUT,  AN  UPPER  TRIANGULAR  ( COMPLEX)  MATRIX,  ON  OUTPUT, 

THE  EXPONENTIAL  OF  B,  WITH  THE  ORIGINAL  B  STORED  IN  THE 
LOWER  TRIANGULAR  PART  (DIAGONALS  ARE  STORED  IN  W(  1 )  TOW(N)) 

W  COMPLEX  ARRAY,  WORK  SPACE 

OVFT  COMPUTER  (REAL)  OVERFLOW  THRESHOLD 

I  ERR  OUTPUT  (INTEGER)  ERROR  STATUS  IF  SOME  EIGENVALUE  OF  B  IS  TOO 
LARGE  TO  EXPONENTIATE,  THE  PROGRAM  QUITS  AND  RETURNS  I ERR=- 2 
N  DIMENSION  OF  MATRIX  B 

NB  LEADING  DIMENSION  OF  ARRAY  B 

INTERNAL  VARIABLES 

INTEGER  I , J ,K, L ,M.N4 , N2 ,KSQ  * 

COMPLEX  X.Y, SHI  FT, CDOTU 

REAL  ZERO, DL I M, SCALE , Z I MAX , Z IMIN , AO , A 1  , A2 , D I  ST 
REAL  DMAX ,  TEMP  ,  EM  I N ,  PMAX  ,  BMAX  ,  AB  ,  G 

DATA  AO, A1  ,A2  /-  1  . 797  804884  ,  1  9  5  8  4  1  4  6  4  0,-0293902  43  90/ 

DATA  DL IM,  ZERO/ 5 , 0 , 0  0/ 

AB(X)=ABS(REAL(X) )+ABS( IMAG(X) ) 

G(  I  )=MAX(  ZERO,  AO-t- 1  •  ( A1  +  I  •  A2  )  ) 

IF  REAL ( 8 ( N , N ) )  IS  TOO  LARGE,  RETURN  I  ERR=- 2  AND  TERMINATE 

I F ( REAL ( B ( N , N ) )  GT  LOG(OVFT) )  THEN 
I  ERR=-  2 
RETURN 
END  I  F 

IF  N=i  OR  N=2 ,  USE  SPECIAL  FORMULA 

I F ( N  LT  1)  RETURN 
I  F ( N  EQ  1 )  THEN 
W(  1 )=B( 1,1) 

B(  1  .  I  )=EXP(W(  1  )  ) 

RETURN 
END  IF 

IF(N. EQ  2 )  THEN 
W(  1 )=B( 1,1) 

W(  2 )=B( 2,2) 

B(  1  ,  1  )=EXP(W(  l )  ) 

B( 2 ,  2 )=EXP(W(  2 )  ) 

B(  2,  1  )=B(  1  , 2 ) 


OOOO  —  O  O  O  -  ooo>  000*i  OOOO 
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I  F  (W(  1  )  EQ  VV(  2  )  )  THEN 
B (  1  , 2 )=B(  I  . 2 ) • B (  1  ,  1  ) 

ELSE 

X=CMPLX(0  5  )  •  (W(  1  )+W(  2  )  ) 

Y=CMPLX( 0 . 0 , 0  5  )  •  (W(  2  )  -W(  1  ) ) 

B ( 1  . 2 )=B( 1  ,  2  )  »  EXP (X)»(S1N(Y)/Y) 

END  IF 
RETURN 
END  IF 

COMPUTE  AND  STORE  THE  EXPONENTIAL  OF  B  IN  THE  LOWER  TRIANGULAR  B  THE 

D I  AGONAL  OF  B  l S  STORED  INW(l),  W( 2 )  ,  ,W(N|  . 

N2=N*N 
N4=N2+N2 
Z  IMAX=0  0 
Z  IMIN=0  0 
DO  2  0  1=1  ,  N 

W(  I  )=B  ( 1,1) 

TEMP=IMAG(W(  I  )  ) 

Z  I  MAX=MAX  (  Z  I  MAX  ,  T  EMP  ) 

ZIMIN=MIN(  ZIMIN,  TEMP-) 

0  CONTINUE 

CALL  CCOPY ( N ,  W,  1  ,  W(  N4+  1  )  ,  1  ) 

SHIFT  OR  SCALE  DCWN  B  TO  GUARANTEE  THE  IMAGINARY  PARTS  ARE  BOUNDED  BY  DLM 

SH I FT=CMPLX( 0.0,00) 

SCALE=1 

KSQ=0 

IF(MAX( -ZIMIN, ZIMAX)  GT.DLIM)  THEN 

I F ( Z IMAX . LT  ZERO  OR . Z IMIN  GT  ZERO )  SH I FT= 

X  CMPLX( 0  0 , (ZIMAX+ZlMIN)/2) 

DO  6  0  1=1  , N 

0  .  W(  I  )=W(  I  ) -SHIFT 

D I  ST=Z  IMAX  -  ZIMIN 
0  I F(DI ST  GT.DLIM)  THEN 

SCALE=SCALE/2  0 
D I ST=D I  ST/ 2  0 
KSQ=KSQ+ 1 
GOTO  80 
END  IF 

IF(KSQ.GT.O)  THEN 
DO  12  0  J=1  ,  N 

CALL  CSSCAL ,( J -  1  , SCALE , B( 1  ,  J  ),  1  ) 

20  CONTINUE 

CALL  CSSCAL ( N , SCAL E , W,  1 ) 

END  IF 
END  IF 

COMPUTE  THE  EXP  OF  B  FROM  BOTTOM  TO  TOP  (RESULTS  STORE  IN  THE  LOWER  B) 

DO  320  I=N, 1,-1 

DO  M0  J=N,  1,-1 

DIST=AB(W(  l  )  -W(  J )  ) 

I F ( D I  ST  LE  G( J -  I  ) )  GOTO  160 
10  CONTINUE 

COMPUTE  THE  SCALE  DIVIDED  DIFFERENCES  AT  W( I )  ,W(N),  RESULTS  STORE 

IN  W(  N+ I  )  ,  ,W(  2N) 

160  CALL  CDDEXP (W(  1  )  ,  W{  N+  I  )  ,  W{  N2+ 1 )  ,  W(  N2+N+ 1  )  .  J -  1+ 1  , OVFT ) 

DO  18  0  K=J+  1  ,  N 

180  W(  K+N )=CMPLX( REAL ( K-  I  )  ) • (W(K+N) -W(K+N- 1 )  ) / (W(K) -W(  I  ) ) 


OOOw  OOOOO^M  OOO  MM  M  OOOM 
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DO  200  K=1 , J 

200  W(K)=W(N4+K)  -  SHIFT 

CALL  CSSCAL (  J  -  1+ 1  , SCALE ,W(  1 )  , 1 ) 

B(  I  .  I )=EXP(  W(  1  ) ) 

I F ( I . EQ.N)  GOTO  320 

DMAX=0  .  0 

DO  2  20  K=N+ 1  ,N2 

2  0  DMAX=MAX(DMAX,AB(W(K)  )  ) 

COMPUTE  E(  1,1+1) . E(l.N)  USING  NEWTON  POLYNOMIAL  METHOD 

Y=W(  N+ 1+  1  ) 

DO  240  10=1+1  ,  N 
W(  K+N2 )=B(  I  ,K) 

B{K,  I )=Y *W(  K+N2 ) 

40  CONTINUE 

DO  3  00  M=I+  2  ,  N 
Y=W(M*-N) 

EMI  N=OVFT 
BMAX=AB(W(M-  1  )  ) 

FMAX=0 . 0 
DO  280  K=N,M, -  1 

X=W(K+N2  )  •  (W(K)  -W(M-  1  )  )  + 

X  CDOTU ( K-Mf  1  ,  W(  N2+M-  1  )  , 1  , B ( M-  1  , K )  . 1 ) 

TEMP=M-  I 

X=X/TEMP 

W(  K+N2  )=X 

B(K,  1  )=B(K,  I  )+X»  Y 

PMAX=MAX  (  PMAX  ,  AB  ( X )  ) 

EMIN=MIN(EMIN.AB(B(K,  I  )  )  ) 
BMAX=MAX(BMAX,AB(W(K)  )  ) 

DO  2  60  L=M-  1  ,K-  1 

60  BMAX=MAX  (  BMAX  ,  AB  (  B  (  L  ,  K )  )  ) 

80  CONTINUE 

TEMP=-  1  0 

I  F  (  BMAX  LT  .  1  .  0  )  TEMP=EMIN+PMAX*DMAX*BMAX 

TEST  FOR  AN  EARLY  CONVERGENCE 

IFfTEMP  EQ.HMIN)  GOTO  320 
!F(PMAX.EQ  ZERO)  GOTO  320 
0  0  CONT I NUE 

20  CONTINUE 

.  END  OF  COMPUTING  THE  EXPONENTIAL  OF  THE  SCALED  B  •••• 

TRANSPOSE  B  AND  RESTORE  THE  ORIGINAL  B  IN  THE  LOWER  PART 

SCAL  E=1  . 0/ SCALE 
DO  330  1=1  ,N-  1 

DO  3  30  J=I+1  , N 
X=B (  I  , J  ) 

B(  I  ,  J )=B ( J  ,  I  ) 

B (  J .  I  )=X* SCALE 
30  CONTINUE 

SQUARE  B  KSQ  TIMES 

DO  420  K=1 , KSQ 

DO  3  4  0  J=1  ,  N 
W(  J  )=W(  J  )+W(  J  ) 

DO  4  0  0  J=N ,1,-1 

CALL  CCOP Y (J,B(l,J),l,W( N+ 1 )  , 1 ) 

CALL  CSCAL ( J ■ 1  ,  W(  N+  J )  , B( 1  ,  J )  , 1 ) 


340 


O  O  O 
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DO  3  80  1=1  ,  J -  1 

CALL  CAXPY (  1  , W( N+ 1 )  , B ( 1 
380  CONTINUE 

B ( J ,  J )=EXP(W{  J ) ) 

400  CONTINUE 

20  CONTINUE 


l).l.B(l.J).l) 


SHIFT  BACK  EXP(B)  AND  RESTORE  THE  DIAGONAL  OF  B  IN  \V(  l  )  TO  W< 


X=EXP (SHI  FT ) 

CALL  CCOPY ( N , W( N4-f- 1  )  ,  1  ,W,  I  ) 
DO  4  80  J=1  ,N 

CALL  CSCAL ( J,X,B(l,J).l) 
460  CONTINUE 
RETURN 
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CDDEXP 

(Subprogram  of  CEXPNT) 


1.  PURPOSE 

The  Fortran  77  subroutine  CDDEXP  computes  the  scaled  divided  differences  at  data 
z(l),  z(2),...,z(n).  They  are  used  in  computing  the  matrix  exponential  by  Newton  poly¬ 
nomial  method.  See  CEXPNT. 


2.  USAGE 


(A).  Calling  Sequence. 

SUBROUTINE  CDDEXP(z,d,r,s,n,ovft) 

Parameters: 

z  is  a  complex  one-dimensional  input  variable  with  dimension  n.  It  con¬ 

tains  the  complex  data  at  which  the  scaled  divided  differences  are  com¬ 
puted. 

d  is  a  complex  one-dimensional  output  variable  with  dimension  n.  On  out¬ 
put,  d  contains  the  scaled  divided  differences  at  z(I),  z(2),...,  z(n). 

r,s  are  two  complex  one-dimensional  variables  of  order  n.  Work  space. 

n  is  an  input  integer  equal  to  the  dimension  of  arrays  z,d,r  and  s. 


ovft  is  the  biggest  positive  real  number,  overflow  threshold. 


S.  DISCUSSION  OF  METHOD  AND  AL  GORITHM 

The  scaled  exponential  divided  differences  at  z(l),...,z(n)  (see  [l])  is  the  first  row  of  the 
exponential  of  a  bi-diagonal  matrix  G,  ={  G,(i,i)=z(i),  i=l,...,n,  G,(i,i+  1)— i, 

i=l . n-l}.  Here  we  apply  a  variation  of  the  scaling  and  squaring  method  (See  [2])  to 

compute  exp(G*).  We  first  scale  down  the  diameter  of  (z(i)}  (by  multipling  2'*)  to 
about  0.7.  Then  we  use  the  Taylor  serie  to  compute  exp(Gz).  Finally  we  scale  and 
square  exp(G7)  k  time.  Because  the  the  special  struction  of  G,,  our  program  keeps 
down  storage  needs  to  a  minimum. 

For  high  accuracy,  the  data  z(l),...,z(n)  must  be  ordered  by  their  real  parts:  Re(z(l)) 
<=*  Re(z(2))  <=  ...  <=  Re(z(n)). 

There  are  approximately  75  executable  Fortran  statements. 


4.  REFERENCES 


[l]  A.  McCurdy,  K.C.  Ng  and  B.N.  Parlett,  Accurate  Computation  of  Divided 
Differences  of  The  Exponential  Function,  Mathematics  of  Computation, 
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)ute  The  Exponential  of 
801-836. 
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5.  PROGRAM  LISTING. 


SUBROUT  1 NE  CDDEXP ( Z , D , S , R , N , OVFT ) 

COMPLEX  Z(N) ,D(N) ,S(N) ,R(N) 

REAL  OVFT 
INTEGER  N 

THIS  SUBROUTINE  COMPUTES  THE  (COMPLEX)  SCALED  EXPONENTIAL  DIVIDED 
DIFFERENCES  AT  Z( 1 ),..., Z(N)  BY  SCALING  AND  SQUARING. 

WE  ASSUME  RE( Z( 1 ) )  <=  ...  <=  RE(Z(N) ) . 

CODE  IN  F77  BY  KWOK-CHOI  NG ,  UC  BERKELEY ;  REV  I S ED  ON  S/18/85 

F7  7  GENERIC  FUNCTIONS  :  ABS  ,  EXP,  LOG,  MAX,  REAL,  IMAG 
F77  NON-GENERIC  FUNCTION:  CMPLX 

REQUIRED  BLAS  SUBROUTINES.  CCOPY ,  CSSCAL ,  CSCAL 
GLOSSARY. 

Z  INPUT  COMPLEX  ABSCISSAE 

D  OUTPUT  COMPLEX  SCALED  EXPONENTIAL  DIVIDED  DIFFERENCES 

S  ,  R  COMPLEX  WORKING  SPACE 

OVFT  COMPUTER  (REAL)  OVERFLOW  THRESHOLD 

N  DIMENSION  OF  ARRAY  - Z , D , S , R 

INTERNAL  VARIABLES 

COMPLEX  SH I  FT , X , Y , C 
INTEGER  I.J.K 

REAL  T , T 1 , T2 , TE , SCALE , RAD  I U , RL IM 
DATA  RLIM/0  7E0/ 

SHIFT  AND  SCALE  DOWN  THE  ABSCISSAE  AND  DETERMINE  THE  NUMBER  OF  SQUARING 

SHIFT  =  CMPLX( 0.0) 

DO  10  1=1  ,N 
0  SHI FT=SH I FT+Z ( 1 ) 

T=N 

SHI  FT=SH  I  FT/T 

IF  EXP (Z(N)-SHIFT)  OVERFLOWS,  THEN  SET  RE( SH I FT)=*E( Z (N) ) . 

IF(REAL( Z(N) -SHIFT) . GT . LOG( OVFT) )  SHIFT= 

XCMPLXf  REAL ( Z (N) ) ,  IMAG( SHIFT) ) 

DO  20  1=1 ,N 
0  Z( I  )=Z( I ) - SH I  FT 
RAD  I  U=0 
DO  3  0  [=1  ,  N 

0  RAD  I U=MAX ( RAD  I U , ABS ( Z ( I ) ) ) 

IF  (RADIU  LE  RL IM)  THEN 
K=0 

ELSE 

K=l+ ( LOG ( RAD  I U )  +  0  3566749)/0  .  6931472 
END  IF 

SCALE  DOiVN  Z.  Z  :  =  Z/2»*K  SUCH  THAT  |  Z  |  <0  7  ,  INITIALIZE  S  AND  D 


T=2  •  •  K 

SCALE=1  0 /T 

RAD IU=«ADI UPSCALE 

40 

DO  4  0  1=1  ,N 

D (  I  )=CMPLX(  1 

0) 

so 

DO  SO  1=1  ,N 

R (  I  )=CMPLX(  0 

0) 
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CALL  CCOPY ( N , D , 1 ,S, 1 ) 

CALL  CSSCAL(N, SCALE, Z , 1 ) 

C 

C  SUMMING  THE  TAYLOR  SERIES  ;  WHILE  J=0  ,  1  ,  .  .  UNTIL  CONVERGE  DO 
C 

J=0 

Tl=EXP( -RAD1U) 

TE=1 

6  0  J=J+ 1 

T=J 

TE=TE*RAD I U/T 
T2=Tl+TE 

I F ( Tl . NE  T2 )  THEN 
X=Z  ( 1 ) • S ( 1 ) 

S( 1 )=X/T 
DO  70  1=2  ,N 
T=I  -  1 

C=Z( 1  )  •  S ( I  )  +  S ( I  -  1  )  *T 
T=J+T 
S(  I  )=C/T 
D (  I  )=©(  I  )+S(  I  ) 

70  CONTINUE 

GOTO  60 
END  IF 

D ( I )=£XP(Z( 1 ) ) 

C 

C  NOW  SQUARE 
C 

100  CONTINUE 

IF(K.EQ  0)  THEN 
X=EXP  (  SH  I  FT  ) 

DO  110  1=1  ,N 

110  Z (  I  )=Z(  I  l  +  SHIFT 

CALL  CSCAL ( N , X , D , 1 ) 

ELSE 

SCALE=1  0 
S  (  1  )=0(  2  ) 

DO  140  1=2  , N 

X=S(  1  ) 

S( 1 )=D( I ) 

S( I )=EXP(Z( I ) ) 

R  (  1  )=(D(  1  )  +  S  (  I  )  )  ♦  D(  I  ) 

DO  130  J=2  ,1-1 
Y=S(  J) 

C=S  ( J -  1  )  •  ( Z (  I  ) - Z ( J  -  1  ) ) 

T=I  -  1 
C=C+X«T 
T=J  -  1 
S( J )=C/T 

R ( I  )=R( I )+S( J ) • D( J ) 

X=Y 


130 

CONTINUE 

SCAL  E=SCALE/  2  0 

R ( 1 )=ft( I ) ‘SCALE 

140 

CONTINUE 

- 

DO  ISO  1=1 ,N 

1  &0 

z (  1  )=Z(  I  )+Z (  I  ) 

CALL  CCOPY ( N , R , 1 , D, 1 ) 
D( I )=EXP(Z( 1 ) ) 

K=K-  1 
GOTO  100 
END  IF 
RETURN 
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CFMUL  V 


1.  PURPOSE. 

The  Fortran  77  subroutine  CFMULV  computes  the  matrix  product  from  right  to  left: 

PFP*V. 

Here  F  is  complex  triangular. 


2.  USAGE. 

(A).  Calling  Sequence. 

SUBROUTINE  CFMUL VJF.P.V.w.n.of.np.nv.k) 

Parameters  : 

F  is  an  input  complex  two-dimensional  variable  with  row  dimension  nf 
and  column  dimension  at  least  n.  F  contains  a  complex  triangular 
matrix  of  order  n. 

P  is  an  input  complex  two-dimensional  variable  with  row  dimension  np 
and  column  dimension  at  least  n.  P  contains  an  unitary  matrix. 

V  is  a  complex  two-dimensional  variable  with  row  dimension  nv  and 
column  dimension  at  least  k.  On  output,  V  is  overwritten  by  the  pro¬ 
duct  PFP*V. 

w  is  a  complex  one-dimensional  variable  with  dimension  at  least  n.  Work 
space. 

n  is  a  input  integer  equal  to  the  order  of  the  matrix  S. 

nf,np  are  input  integers  equal  to  the  row  dimension  of  F  and  P  respectively. 

nv,k  are  input  integers  equal  to  the  row  and  column  dimension  of  V  respec¬ 
tively. 


S.  DISCUSSION  OF  METHOD  AND  ALGORITHM. 

There  are  approximately  12  executable  Fortran  statements. 


4.  REFERENCES. 


NONE 


ooo*  «  ooo  ooo  ooooooooooooooooooo 


SUBROUTINE  CFMULV { F , P , V , W,  N , NF , NP , NV , K) 

COMPLEX  F(NF,N) ,P(NP,N) ,V(NV,K) ,W(N) 

INTEGER  N , NF , NP , NV , K 

THIS  SUBROUTINE  EVALUATES  V  : =  P • ( F • ( PT* V ) ) 

CODE  IN  F77  BY  K.C  NG ,  UC  BERKELEY;  REVISED  ON  6/8/85. 

REQUIRED  BLAS  SUBROUTINE  :  CAXPY 
REQUIRED  BLAS  FUNCTION  CDOTC 

GLOSSARY 

F  COMPLEX  TRIANGULAR  MATRIX 

P  UNITARY  MATRIX 

V  N  BY  K  MATRIX:  ON  OUTPUT,  IT  IS  REPLACED  BY  P»(F»(PT«V)) 

W  COMPLEX  ARRAY,  WORK  SPACE 

N  DIMENSION  OF  MATRICES  F,P,V 

NF  LEADING  DIMENSION  OF  ARRAY  F 

NP  LEADING  DIMENSION  OF  ARRAY  P 

NV  LEADING  DIMENSION  OF  ARRAY  V 

INTERNAL  VARIABLES: 

INTEGER  I , J 

COMPLEX  CDOTC,  T,  ZERO 
DATA  ZERO / ( 0E0 , 0E0 )  / 

DO  ONE  COLUMN  AT  A  TIME 

DO  100  J—l , K 

FORM  W  =  F  •  (  PT  •  V  (  .  ,  J  )  ) 

DO  2  0  1=1  ,  N 

0  W(  I  )=ZERO 

DO  <0  1=1  ,  N 

T=CDOTC ( N , P ( 1 , I ) , 1 , V{ 1 , J ) , 1 ) 

CALL  CAXPY (  I  , T , F (  1  ,  I  )  ,  1  ,  W,  1  ) 

0  CONTINUE 

FORM  V ( . , J )  =  P  •  W 


60 


DO  60  1=1 , N 
V (  I  ,  J  )=ZERO 
DO  8  0  1=1  ,N 


bias  c 


BLASJC 


1.  PURPOSE. 

The  Basic  Linear  Algebra  Subroutines  (Complex)  compute  the  basic  matrix-vector 
operation.  See  [l]. 


2.  USAGE.  .. 

(A).  Calling  Sequence. 

SUBROUTINE  CAXPY  (N,CA,CX,INCX,CY,INCY) 
SUBROUTINE  CCOPY  (N,CX,INCX,CY,INCY) 
SUBROUTINE  CSSCAL(N,SA,CX,INCX) 

SUBROUTINE  CSCAL  (N,CA,CX,INCX) 

COMPLEX  FUNCTION  CDOTU  (N,CX,INCX,CY,INCY) 
COMPLEX  FUNCTION  CDOTC  -  (N,CX,INCX,CY,INCY) 

Parameters  : 

See  listing. 


3.  DISCUSSION  OF  METHOD  AND  ALGORITHM. 

See  [  1] . 

There  are  approximately  75  executable  Fortran  statements. 


4-  REFERENCES. 

[lj.  J.J.  Dongarra,  C.B.  Moler,  J.R.  Bunch,  and  G.W.  Stewart,  UNPACK  users' 
guide,  SLAM,  Philadelphia/1979. 


oooo  o  oooo  ooo  oooo  o  oooo  ooooo 


5.  PROGRAM  LISTING. 


UNPACK  BLAS  FUNCTIONS  AND  SUBROUTINES 
CDOTC ,  CDOTU 

CAXPY ,  CCOPY ,  CSSCAL,  CSCAL 


SUBROUTINE  CAXPY ( N , CA , CX , INCX.CY, INCY) 

CONSTANT  TIMES  A  VECTOR  PLUS  A  VECTOR. 

JACK  DONGARRA.  LINPACK,  3/11/78 

COMPLEX  CX( 1 ) ,CY( 1 ) ,CA 
INTEGER  I , INCX, INCY, IX, IY,N 

I F ( N . LE . 0 ) RETURN 

IF  ( ABS ( REAL ( CA ) )  +  ABS ( AIMAG ( CA ) )  EQ.  00  )  RETURN 
I F ( I NCX  EQ  1  AND  I NCY  EQ .  1  )  GO  TO  20 

CODE  FOR  UNEQUAL  INCREMENTS  OR  EQUAL  INCREMENTS 
NOT  EQUAL  TO  1 


IX  =  1 
IY  =  1 

I  F  (  I  NCX .  LT  0  )  I X  =  (  -  N+- 1  )  *  I  NCX  +  1 
IF(  INCY.  LT  0)IY  =  (-N+-1  )•  INCY  +  1 
DO  10  I  =  1  ,N 

CY(IY)  =  CY( I Y)  +  CA*CX( IX) 

IX  »  !X  +  INCX 
IY  -  IY  +  INCY 
10  CONTINUE 
RETURN 

CODE  FOR  BOTH  INCREMENTS  EQUAL  TO  1 

20  DO  30  I  =  1 ,N 

CY ( I )  =  CY( I )  +  CA*CX( I ) 

30  CONTINUE 
RETURN 
END 


SUBROUTINE  CCOPY{ N , CX , I NCX , CY , INCY) 

COPIES  A  VECTOR,  X ,  TO  A  VECTOR,  Y. 

JACK  DONGARRA,  LINPACK,  3/11/78 

COMPLEX  CX( 1 ) ,CY( 1 ) 

I  NT EGER  I  ,  I NCX ,  INCY ,  I X ,  I Y , N 

I F ( N  LEO) RETURN 

I F ( INCX.EQ. 1  AND  INCY  EQ  1 )GO  TO  20 

CODE  FOR  UNEQUAL  INCREMENTS  OR  EQUAL  INCREMENTS 
NOT  EQUAL  TO  1 

IX  *  1 
IY  -  1 

IF( INCX. LT. 0 ) IX  =  ( - N+ 1 ) • I NCX  +  1 
I F( INCY . LT. 0 ) IY  =  ( -N+l ) • INCY  +  1 
DO  10  I  »  1  ,N 


ooooo  non  noon  o  oooo  ooo 
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CY( IY)  =  CX( IX) 

IX  =  IX  +  INCX 
IY  =  IY  +  INCY 
10  CONTINUE 
RETURN 

CODE  FOR  BOTH  INCREMENTS  EQUAL  TO  1 

20  DO  30  I  =  1 , N 

CY (  I  )  =  CX(  I  ) 

30  CONTINUE 
RETURN 
END 


COMPLEX  FUNCTION  CDOTU ( N , CX , INCX , CY , INCY ) 

FORMS  THE  DOT  PRODUCT  OF  TWO  VECTORS . 

JACK  DONGARRA,  L INPACK,  3/11/78. 

COMPLEX  CX ( 1 ) , CY ( 1 ) , CTEMP 
INTEGER  I , INCX, INCY, IX, IY,N 

CTEMP  =  (00,00) 

CDOTU  =  (00,0.0) 

I F ( N . LE  0 ) RETURN 

I F ( INCX.EQ. 1  AND  INCY.EQ. 1 )GO  TO  20 

CODE  FOR  UNEQUAL  INCREMENTS  OR  EQUAL  INCREMENTS 
NOT  EQUAL  TO  1 

IX  =  1 
IY  =  1 

I F ( INCX. LT. 0) IX  =  ( -Nf  1 )* INCX  +  1 
I F ( INCY  LT  0) IY  =  (-N*l )*INCY  +  1 
DO  1 0  I  =  1  ,  N 

CTEMP  =  CTEMP  +  CX( IX) *CY( I Y) 

IX  =  IX  +  INCX 
IY  =  IY  +  INCY 
10  CONTINUE 

CDOTU  =  CTEMP 
RETURN 

CODE  FOR  BOTH  INCREMENTS  EQUAL  TO  1 

20  DO  30  I  =  1 ,N 

CTEMP  =  CTEMP  +  CX  (  I  )  •  CY (  I  ) 

30  CONTINUE 

CDOTU  =  CTEMP 

RETURN 

END 


COMPLEX  FUNCTION  CDOTC ( N, CX , I NCX , CY , INCY ) 

FORMS  THE  DOT  PRODUCT  OF  TWO  VECTORS,  CONJUGATING  THE  FIRST 
VECTOR . 

JACK  DONGARRA,  L INPACK,  3/11/78. 

* 

COMPLEX  CX( 1 ) ,CY( 1 ), CTEMP 
INTEGER  I , INCX, INCY, IX, I Y,N 

CTEMP  =  (0.0,0  0 ) 

CDOTC  =  (0.0,0  0 ) 
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IF (N. LE . 0 ) RETURN 

IF( INCX. EQ. 1  AND  INCY. EQ  1 )G0  TO  20 

CODE  FOR  UNEQUAL  INCREMENTS  OR  EQUAL  INCREMENTS 
NOT  EQUAL  TO  1 

IX  =  1 
IY  =  1 

IF(  INCX.LT.  0  )  IX  =  (  - N-t- 1  )  •  I NCX  +  1 
IF(  INCY  LT.  0  )  IY  =  ( -N+-1  )•  INCY  +  1 
DO  10  I  =  1  ,  N 

CTEMP  =  CTEMP  +  CONJG(CX( IX) ) *CY( IY) 

IX  =  IX  +  INCX 
IY  =  IY  +  INCY 
10  CONTINUE 

CDOTC  =  CTEMP 
RETURN 

CODE  FOR  BOTH  INCREMENTS  EQUAL  TO  1 
20  DO  30  I  =  1  ,  N 

CTEMP  =  CTEMP  +  CONJ G(CX( I ) ) *CY( I ) 

3  0  CONTINUE 

CDOTC  =  CTEMP 

RETURN 

END 


SUBROUTINE  CSSCAL ( N , SA , CX , INCX ) 

SCALES  A  COMPLEX  VECTOR  BY  A  REAL  CONSTANT 
JACK  DONGARRA,  UNPACK,  3/U/78. 

MODIFIED  FOR  F77  BY  K.C.  NG ,  4/H/85. 

COMPL EX  CX( 1 ) 

REAL  SA 

INTEGER  I , INCX.N.NINCX 

IF(N  LE  0 ) RETURN 
I F ( INCX  EQ. 1 )GO  TO  20 

CODE  FOR  INCREMENT  NOT  EQUAL  TO  1 

NINCX  =  N* INCX 
DO  10  1=1 .NINCX, INCX 
CX( I )  =  CX( I ) • SA 
RETURN 

CODE  FOR  INCREMENT  EQUAL  TO  1 

DO  3  0  I  =  1  ,  N 
CX ( I )  =  CX( I ) • SA 
RETURN 
END 


SUBROUTINE  CSCAL ( N , CA , CX . I NCX ) 

SCALES  A  VECTOR  BY  A  CONSTANT 
JACK  DONGARRA,  L INPACK,  3/11/78 

COMPLEX  CA , CX ( 1 ) 

INTEGER  I , INCX.N.NINCX 


IF (N. LEO ) RETURN 
I F ( INCX.EQ. 1 )GO  TO  20 

CODE  FOR  INCREMENT  NOT  EQUAL  TO  1 

NINCX  =  N»  I NCX 
DO  10  1=1 .NINCX, I NCX 
CX( I )  =  CA*CX( I ) 

10  CONTINUE 
RETURN 

CODE  FOR  INCREMENT  EQUAL  TO  1 

2  0  DO  3  0  l  =  1  , N 

CX( I )  =  CA*CX( I ) 

30  CONTINUE 
RETURN 


Legal  Notice 


This  report  was  prepared  as  an  account  of  work  sponsored  by 
the  Center  for  Pure  and  Applied  Mathematics.  Neither  the 
Center  nor  the  Department  of  Mathematics,  makes  any  warranty 
expressed  or  implied,  or  assumes  any  legal  liability  or 
respons ibi 1 i ty  for  the  accuracy,  completeness  or  usefulness 
of  any  information  or  process  disclosed. 


